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Abstract. The growth of solid particles towards meter sizes in protoplanetary disks has to circumvent at least 
two hurdles, namely the rapid loss of material due to radial drift and particle fragmentation due to destructive 
collisions. In this paper, we present the results of numerical simulations with more and more realistic physics 
involved. Step by step, we include various effects, such as particle growth, radial/vertical particle motion and 
dust particle fragmentation in our simulations. We demonstrate that the initial dust-to-gas ratio is essential for 
the particles to overcome the radial drift barrier. If this value is increased by a factor of 2 compared with the 
canonical value for the interstellar medium, km-sized bodies can form in the inner disk (< 2 AU) within 10^ yrs. 
However, we find that solid particles get destroyed through coUisional fragmentation. Only with the unrealistically 
high-threshold velocities needed for fragmentation to occur (> 30 m/s), particles are able to grow to larger sizes 
in disks with low a values. We also find that less than 5% of the small dust grains remain in the disk after 1 Myrs 
due to radial drift, no matter whether fragmentation is included in the simulations or not. In this paper, we 
also present considerable improvements to existing algorithms for dust-particle coagulation, which speed up the 
coagulation scheme by a factor of ~ 10*. 
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1. Introduction 

The coagulation of sub-/im dust particles is believed 
to be the initial step of planetesi mal formation in 



of planetesima ls by particle coagulation (jYoudinl 12004 : 
Dominik et ah 2007; Brauer ct al. 20"07[). 

One of these obstacles is the radial inward drift of 
solid bod i es tow ards the central star as first dicussed by 
Whipple! (|l972l ) and IWeidenschillind (|l977al ). The gas in 



disks around pre-inain-se quence stars (jKlahr fc Brandner 



20061 iNatta et al.l |2007|) . Evidence for 
that are found 



dust grain evo- 
n the interstel- 



lution beyond sizes ttiat are 

lar medium is provided by mid-in frared spectroscopy of 



disks around Herbig Ae/Be stars (IBouwman et all 200ll; 



van Boekel et al. 20 031), T Taur i stars (jPrzygodda et al.l 



2003t iKessler-Silacci et al.l 12007) and also around brown 
dwar fs (JApai et al.ll2004 l2005l l2007t ISiciha-Aguilar et al ' 
20071) . Millimeter interferometry indicates 
lations of particles which have grown to 



large 
even 



Wilner et al.l 120051 iRodmann et al.ll2006[ ) 



popu- 
larger 



sizes, ranging up to several centimeters ( Test! et al.ll2003t 



All these observations give reason to model the evo- 
lution of particles in protostellar disks in order to ex- 
plain t he observations,! data (iDullcmond & Dominik 2004, 
20051 iTanaka et all l2005t [Nomura fc Nakagawal l2006t 



D'Alessio et al.l l2006t lOrmel fc Cuzzil l2007t) . However, 
these theoretical investigations do not only attempt to 
model the evolution of the appearance of protostellar 
disks. They also unveil certain obstacles in the formation 



a protostellar disk moves slightly sub-Keplerian due to a 
radial pressure gradient. For this reason, the dust which 
moves with near Keplerian velocity feels a continuous head 
wind of gas. Hence, the particle loses angular momentum 
due to drag forces between gas and dust and spirals in- 
ward. If this radial drift of the particles is not prevented 
by some mechanism, then the solid particles drift into the 
inner evaporation zone and are lost for the process of plan- 
etesimal formation. To give an example, the radial drift 
time scale for meter-sized bodies at 1 AU is ^ 10^ yrs. 
Within this time scale these boulders at 1 AU drift into 
the inner regions of the disk and evaporate. A possible way 
out of this problem is particle growth since the radial drift 
velocity is fairly dependent on the particle radius. For ex- 
ample, the drift velocity of meter-sized particles at 1 AU in 
the disk is ~ 50 m/s, but the radial drift velocity of 10 m 
sized bodies is already 10 times lower. Therefore, swift par- 
ticle growth could prevent the particles from drifting into 
the evaporation zone. However, the general disk evolution 
comprises a considerable particle loss due to evaporation 
which is hard to prevent. This problem is a major topic of 
this paper. Other su blimation zones o f the disk, e.g. the 
snow line at '--^ 2 AU (JLecar et al.ll2006l) , could also play a 
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role for particle drift and coagulation processes. However, 
we will for now neglect this issue which will be the topic 
of a forthcoming paper. 

Another obstacle is the fragmentation of solid parti- 
cles. While low-velocity collisions lead to part icle growth 
high velocity impa,cts lead to destruction (IBlum et al 



1998t IPoppe et al.l Il999t iBlum fc WuT^ I2OOOI ). For ex 
ample, the relative particle velocity of meter-sized bod- 
ies in a protostellar disk can be more than ^ 3 m/s 
('Weidenschilling Il977at iMarkie wicz et al.l Il99ll ). iBenz 
(2000) found that meter-sized rocks appear unlikely to 
survive an impact with a relative low collision velocity of 
some cm/s. For porous objects, collision velocities higher 
than 4% of t he sound speed lead to particle destruction 
( Sirondl2004l ). For this reason, the particle size of roughly 
a meter seems to pose an upper limit for particle coagu- 
lation. 

These two obstacles, the radial drift barrier and the 
fragmentation barrier, are the issue of this paper. We 
present a disk model including the growth, the radial drift 
and the fragmentation of the particles. We show how these 
three effects change the evolution of the disk by including 
them step by step. 

1. In the first step we only consider particle coagulation 
due to Brownian motion, vertical settling and tur- 
bulent mixing. This step shows to which sizes parti- 
cles can grow if radial drift and fragmentation are ne- 
glected. 

2. The second step includes the radial drift and the radial 
mixing of dust. The particles are now allowed to move 
inwards and to disappear into the evaporation zone. 
However, we investigate which disk parameters influ- 
ence the drift time scales and for which parameters the 
dust particles overcome the drift barrier. 

3. The last step also includes particle fragmentation. We 
show under which conditions, i.e. in which regions of 
the disk and for which disk parameters, it is possible 
for the dust to overcome this barrier. 

The radial drift barrier is not only of interest for the radial 
drift itself. Particles close to this barrier are most suscep- 
tible to the motions of the gas and the gravitational effects 
of the dust. For example, particles can be trapped in very 
elongated gas vortices in magnetorotational tur bulence 
(JBalbus fc Hawlevlll99lHBarge fc Sommeriall 19951 ). These 
effects can slow down the radial drift by a factor of two 
( Johansen et al.ll2006l ). Under certain conditions the solid 
particle layer itself may b ecome gravitationally unstable 
(| Johansen fc YoudinI 120071 ). In high dust density regions, 
the particles contract due to their own g ravity and may 



form a planetesimal within a few orbits (jJohansen et al 



20071). Moreover, the flow of the gas and the dust can be 
unsta ble to the streaming instablity (JYoudin fc Goodman 
20051) which leads to particle clumping, and possibly also 
to a gravitational collapse of the dust. Apparently, the ra- 
dial drift barrier is not only connected to the radial motion 
of the dust particles, but involves various other important 
effects as well. For this reason, it is vital to answer the 



question if particles can actually reach the size regime at 
which non-linear effects become of importance. 

In this paper we will implement a 2+1 dimensional 
model. The first dimension is the radial coordinate of the 
disk r, the second one is the height above the midplane 
z and the third coordinate is the mass of the dust par- 
ticles m. The dust may move radially due to radial drift 
and radial mixing. We will numerically solve the continu- 
ity equation for this problem for each particle species. The 
time evolution for the particle size distribution is deter- 
mined by the coagulation equation. We will numerically 
solve this equation as well. In the vertical direction, we 
will always assume that each particl e species is in vertical 
sedimentation/mixing equilibr ium (jDubrulle et al.l Il995 ; 
ICuzzi fc Weidenschillingl 120061 ). Hence, we will not solve 
the time dependent continuity equation in the vertical 
direct ion as done for example by 1 Dullemond fc DominikI 
(|2005l ). Nor do we need to solve the coagulation equa- 
tion at all z explicitly. Instead, we solve the vertically in- 
tegrated coagulation equation, which significantly saves 
computational time (cf. Appendix IB|) . We also formulated 
the coagulation equation in an implicit way (cf. Appendix 
ICj) which saves another factor of ~ 100 of computer sim- 
ulation time. 

2. Model equations 

2.1. Disk model 

We consider a disk of mass Afdisk and an inner and an 
outer disk radius of ri„ = 0.03 AU and rout = 150 AU. We 
adopt a central stellar mass of 0.5 Mq and a disk mass 
of 0.01 Mi, if not otherwise noted. The mass distribution 
of gas and dust inside the disk is given by the surface 
density S. This quantity generally depends on the distance 
to the central star r and the azimuthal angle ip. However, 
since we assume a disk which is axisymmetric the surface 
density E only depends on the radius r. We assume that 
this dependency can be described by a power law 



E = E( 



llAuJ 



(1) 



The power law index of the surface densi ty S is set to be 



0.8 in the course of this paper followin g iKitamura et al 



(I2OO2I ) and [Andrews fc William^ (|2007l ). The surface den- 
sity at 1 AU, which is denoted by Sq, is chosen in a way 
that the condition 



2tt 



T.(r)rdr = Mdisk 



(2) 



holds. With the values mentioned above, the surface den- 
sity of the gas at 1 AU in the disk can be calculated to be 
- 20 g/cm^. 

This surface density distribution of gas and dust, on 
which the dust particle coagulation calculations in this 
paper are based, differs significantly from the disk model 
which is usually referred to as the minimum mass solar 
nebula (MMSN) model (|Weidenschillinslll977bt iHavashi 
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Fig. 1. The surface density distribution of the gas as a 
function of disk radius for the disk model discussed in 
the paper at hand and the MMSN model as discussed in 
Sec. [21 Note, that the surface densities at 1 AU differ by 
more than one order of magnitude. 



19811 ). The fundamental difference can be found in the 
distribution of the mass of gas and dust in the disk. In 
the MMSN disk model, the power law index of the surface 
density is S ~ 1.5, while the disk model in this paper 
adopts 6 = 0.8. Fig. [T] shows the surface densities as a 
function of disk location for both disk models assuming a 
disk mass of 0.01 M^.. The MMSN model implies surface 
densities of ~ 600 g/cm^ at 1 AU in the disk. With the 
disk model in this paper, we yields surface density values 
of ~ 20 g/cm^ which is more than one order of magnitude 
lower. 

The actual distribution of mass in a protoplanetary 
nebula is still a matter of debate. There is evidence from 
meteoritics that the densities in the protosolar nebula in 
the planet forming region have been very high, implying 
disk masses much larger than the MMSN (JDesch et al.l 
2002r) . On the other hand, resolved millimeter dust emis- 
sion maps of protoplanetar y nebula seem to indicat e 
much lower surface densities ( Andrews fc Williamsll2007 ). 
However, millimeter dust observations of disks may not 
trace the radial profile of the gas mass density correctly 
since particle growth to larger sizes is expected to proceed 
more quickly in the inner parts of the disk than in the 
outer parts. For this reason, the dust continuum emission 
becomes flat even though the radial profile of the sur- 
face density might have a steep radial behaviour. Hence, 
analysing dust emission maps assuming a constant dust 
particle size throughout the disk likely leads to power law 
indices S which are systematically shifted towards lower 
values. 

The actual surface density distribution probably lies 
in between these two extreme cases, i.e. the MMSN model 
with ^ = 1.5 and the obse rvational median oi S = 0.5 
( Andrews fc Wilhamd 120071 ). The surface density profile 
adopted in this paper is chosen to be between these two 



extremes. Our goal is to gain insight in dust disk evolution 
in a environment which might well be a likely scenario 
considering the S range under discussion. 

We assume that the gas in the disk is in a steady state, 
even for times as long as 1 Myr. Hence, the gas densi- 
ties in our model do not change in time. We only focus 
on the dust component in the disk which evolves on a 
steady gas background. To unveil the robustness of this 
assumption, we compare the following time scales. Particle 
growth time scales are of the order of '--^ 10^ orbital times 
before fragmentation prevents further particle growth (c.f. 
Sec. 13. 31) . Radial gas accretion velo cities are of the order o f 
~ 1 ... 10 cm/s at 1 AU in the disk (JTakeuchi fc Linll2002l ). 
With regard to 1 AU, this leads to accretion time scales 
of the order of 10^ yrs, which is much lar ger than typical 
particle growth time scales. For example, iTakeuchi fc Lira 
(2002) find that in the first lO**"^ yrs, the gas surface 
density between 1 and 100 AU is hardly affected by vis- 
cous evolution. However, after 10^ yrs, the surface den- 
sity of the gas may change signifi cantly over time scales 
of several Myrs ( Reves-Ruidl2007l ) . This introduces a sys- 
tematic uncertainty in our dust evolution model regarding 
late evolutionary stages of T Tauri disks. 

The temperature T is assumed to be the midplane tem- 
perature of a disk irradiated under an angle of ai^ = 0.05 
around a T Tauri star with a surface temperature of 
T^ = 4000 K and a radius of i?* = 2.5 Rq. If we as- 
sume the disk to be isothermal in the vertical direction 
then the temperature is given by 



rr 1/4 

1 ^ a- 

irr 



-1/2 



r 
R. 



(3) 



With this dependency the temperature at 1 AU is 204 K. 
The evaporation temperature at 0.03 AU is ^ 1400 K. 

2.2. Vertical structure of gas 

We consider a thin disk, which means that z <C r. The 
quantity z denotes the height above the midplane. Under 
this condition the vertical mass density distribution of the 
gas can be described by 



Pg{z,r) 



E(r) 



2ttH 



exp ( -~z^/2H^ ) 



(4) 



In this expression the quantity H denotes the pressure 
scale height of the gas given by iJ — Cg/^k, where 
Cg — ^JkT / fi is the isothermal sound speed and f2k — 
yjGMi,/r^ denotes the Kepler frequency. The quantities 
k and G are the Boltzman constant and the gravitational 
constant, respectively. The mean molecular weight /i is as- 
sumed to be 2.3 TOp (mixture of molecular hydrogen and 
helium) where m-p is the mass of a proton. 

With Eq. dU, the gas mass density in the midplane of 
the disk at 1 AU in our model with 5 = 0.8 is given by 
10~" g/cm^. Adopting the MMSN model with 5 = 1.5 
leads to a gas mass density of 4 x 10"^*^ g/cm^ which is 
more than one order of magnitude higher. 
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2.3. Dust variables 

Before we introduce the vertical structure of the dust we 
define some variables that describe the dust distribution. 
We define pj°* to be the total mass of the dust per cm'^ 
at a certain point in space. To describe the particle mass 
distribution we define a dust density pd{rn), normalised 
such that 



tot 



pd{in) dm 



Now, we introduce the number density n by 



n{m) 



Pd{m) 



(5) 



(6) 



This quantity gives the number of solid particles of a cer- 
tain mass m in a, unity mass interval. Now, we can define 
the integrated number density w and the surface density 
of the dust by 



oj{ra) 



+ OC 



n{m) Az 



(7) 



and 



Sd(?Ti) — muj(m) 



(8) 

To implement these expressions in a computer program 
we have to introduce a mass grid {wk} and a measure 
{Awk}. With the definitions of the number density A^k 
and the dust density pk on the mass grid 

iVk — n(TOk)Amk and pk — Pd("ik)Amk (9) 



equation ([6|) implies 



Pk = mkiVk 



(10) 



The quantity n(mk) in Eq. ^ is an arbitrary value of the 
function n within the interval Arrik around my^. To define 
a surface density on a lattice we also have to introduce a 
vertical space grid {^i} and its measure {Azi}. The surface 
density is then given by 



"^k > ^iVk(zi) Azi. 



I 



2.4. Vertical structure of the dust 



(11) 



In a protostellar disk the scale height of the dust is deter- 
mined by an equilibrium between two processes, namely 
the settling of the dust towards the midplane of the disk 
due to vertical gravity, and the vertical mixing of the dust 
due to turbulent diffusion. The more turbulence there is, 
the harder it is for the dust to form a thin layer since it 
is mixed up and transported back to the higher regions of 
the disk. 

To describe turbulence, w e will use the a-prescription 
of lShakura fc SunvaevI (|l973[) . In this prescription the tur- 
bulent diffusion coefficient of the gas at a certain radius r 
in the disk is parameterized by the scale height of the gas 
H and the isothermal sound speed Cs by 



D„ — ac=H 



(12) 



a determines the amount 
Ob servations suggest 



(JHartmann et al.l llQQSf ). 



The dimensionless parameter 
of turbulence in the disk, 
turbulent a value of 10^'^ 
Numer ical simulations of the magneto rotational insta- 
bility (JBalbus fc Hawlevlll991l) yield turbulence parame- 
ters of the same order of i nagnitude (IBrandenburg et al 



a 



P r 



1995; H awlev et al.l 119951 : iSano et al.l 



20041 ) 



Moreover, 
is a minimal 



Weidenschillinsi ( 1980 ) showed that there 
amount of turbulence in every protostellar disk corre- 
sponding to an a value of about 10~^. 

In addition to the a-value we have to introduce a sec- 
ond dimensionless number in order to describe the vertical 
structure of the dust. This so-called Stokes number is de- 
fined by 

(13) 



btk = Sik Ct ^ 



CsPg 

The variable Ok and ps denote the radius of the dust par- 
ticle of mass k and its material density, respectively. If the 
Stokes number is much smaller than unity, then the dust 
particles are strongly coupled to the gas. In this case, the 
motions of the dust are basically the motions of the gas 
and both components have the same behaviour with re- 
gard to diffusion. If St exceeds unity, then the particles 
decouple from the gas and are hardly influenced by the 
turbulent motions of the gas. The turbulence parameter 
q in Eq. (J13p determines whether turbulent diffusion in 
the disk is realized by big turbulent eddies moving slow 
(q = 1) or by small turbulent eddies moving fast (q = 0). 
Throughout this paper we will assume that q = 1/2 follow- 
ing ICuzzi et all (|200ll) and ISchrapler fc Hennind (|2004l ) 
unless otherwise stated. 

The dimensionless turbulence parameter q is also con- 
nected with the velocity vt of the large turbulent eddies, 

vt = a'^Cs , (14) 

which significantly inffuences the relative turbulent ve- 
locities of the dust (cf. Sec. 12.6. ip and, hence, its co- 
agulation and fragmentation time scales. Various au- 
thors have used quite different values for q during 
the past decades which led to very different turbulent 
eddy velocities and, hence, different rela tive par ti cle ve - 
locities produced by tur bulence. While MorfiU ( 19881) . 
Weidenschilih^ (|l988l ) or I Weidenschilling fc Cu"^ (|l993l ) 
use turbulent gas velocities of acs, which implies q = I, 
more recent publications explici tly der i ved q =1/2 which 



leads to Vt = ^/ac■, (Dubrulle et al.l Il995t ICuzzi et al 



20011 : ICuzzi fc Weidenschillinsll2006l ). If q exceeds 1/2 then 
the time scale of the largest eddy becomes larger than 
an orbital time scale since Teddy '^ ci^~^'^/^k- Turn over 
frequencies smaller than the Kepl er frequency appea r un 



physical to us. Therefore, we foUow lCuzzi et al.l ( 2001 ) and 
adopt q = 1/2. 

With the two dimensionless numbers, a and St, the 
scale heigh t of the dust fek of a certain grain mass Wk 



20061 : iBrauer et al.ll2007l) 



is given by (Dubrulle et al. l ll995l : ICuzzi fc Weidenschilhnd 



a 



min(Stk,l/2)(l + Stk 



(15) 
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Since the dust scale height /ik can not exceed the gas scale 
height H we restrict ft,k to be at most H . With this last 
expression the vertical structure of the dust particles with 
mass TOk is given by 

Pk(r, z) = -p- exp ( -z^/2hl ) . (16) 

In this equation Sk denotes the surface density of the dust 
with mass TTik- 

2.5. Radial motion of the dust 

In this section we will present the equations of radial mo- 
tion for solid particles. We first recapitulate the radial drift 
of individual particles and particle motion due to gas ac- 
cretion. After this we introduce the equations for radial 
mixing due to turbulent diffusion. Finally, we will discuss 
the continuity equation. 

The following equations, which describe the radial drift 
of individual dust particles, are valid as long as the dust- 
to-gas ratio e does not exceed unity. If, however, the dust 
density becomes larger than the gas density then the dust 
starts to have a back-reaction on the motion of the gas 
in a non negligible way. One possible scenario, where the 
dust-to-gas ratio can exceed unity, is when the particles 
settles into a thin midplane layer due to low turbulence 
in the disk. In this dense midplane layer these so-called 
collective effects can become of importance and the radial 
drift equations have to be modi fied in an approp r iate w ay. 
This was discussed in detail bv lNakagawa et al.l (|l986l ). 



2.5.1. Radial drift of individual particles 

We consider two different sources for the radial drift of 
solid particles. The first one is the radial drift of individual 
particles itself. The second source is due to the accretion 
process of the gas. 

The dust particles behave entirely independently and 
the gas is assumed not to be affected by the dust at all. The 
crucial particle characteristics, that determine the drift 
of solid particles, is the Stokes number introduced in the 
last section. In other words, the coupling strength between 
the gas and the dust determines the radial drift. With 
this quantity the radial drift of individual dust particles 
of ma. ss TOk is given by (| Whipple! Il972t IWeidenschilline 
1977ah 

^'dust,k = -— — : — i- ■ (1') 



Stk 



1 

Stk 



The radial drift of the dust is maximal for Stk = 1. For 
this reason, the radial drift barrier can be regarded as a 
region around Stk — 1- The maximum drift velocity «„ in 
the last equation can be calculated according to 



2V,\'^\ 



(18) 



With the numbers mentioned in the last section the last 
expression yields a maximum drift velocity of 45 m/s at 
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1 000.0 
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Fig. 2. The total inward radial velocity of a single dust 
particle as a function of Stokes number St (solid line) as 
discussed in Sec. 12.5.11 The dotted line denotes the accre- 
tion velocity of the gas. The turbulence parameter a is 
10"'^ in this calculation. 



1 AU. Since the power law index of the temperature is 
— 1/2 the maximum drift velocity un is independent of 
the location in the disk. 

The second source for radial velocity of the dust is due 
to the accretion process of the gas. In a gaseous disk angu- 
lar momentum is transported to the outer regions of the 
disk under the action of viscous stress. This transport of 
angular momentum is connected with a radial accretion 
velocity of the gas. This gas velocity also affects the mo- 
tions of the dust. Small particles are strongly coupled to 
the motions of the gas and if the gas moves inwards then 
also the small dust particles move inwards. For larger par- 
ticles, i.e. particles with Stokes number larger than unity, 
the motions of the gas become less and less important 
since larger dust particles decouple from the gas. 



Takeuchi fc Lin! (J2002l ) have calculated the radial ac- 



cretion velocity of the gas due to viscous stress, 



= —3a 



Vv V2 



-6 



(19) 



In the last expression we already adopted the temperature 
dependency T ex r~^'^ according to eq. ([3]). The quantity 
Vk denotes the Kepler velocity Vk — ^k?"- With a turbu- 
lence parameter of a = 10""^ and the values that were 
already mentioned in the last section the last equation 
yields a radial gas accretion velocity of 5 cm/s. 

Now, the total ra dial drift ve locity of solid particles of 
mass mk is given by (JKornet et al. 2001 ) 



-„tot _ „ I "gas 

^dust,k — ^dust,k + 



1 + Stk" 

A plot of this total drift velocity is shown is Fig. ^ 



(20) 



6 
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2.5.2. Radial mixing of the dust and the continuity 
equation 

Turbulence can be interpreted as a kind of diffusion. The 
turbulent diflFusion coefficient for the gas Dg was already 
introduced in the last sections. The equivalent quantity for 
the dust of a certain grain mass is given by (IVolk et al.l 
19801: ICuzzi et al.lll993l : ISchrapler fc Hennindbooj) 



^d,k = 



D^ 



1 + Stk 



(21) 



This relation was already implicitly used in the expression 
for the dust layer thickness Eq. (|15p . For small particles, 
i.e. particles with Stk < 1, the diffusion coefficient for 
the dust matches the diffusion coefficient for the gas. For 
Stokes numbers larger than unity Z^d.k decreases contin- 
uously since the particles more and more decouple from 
the turbulent gas. 

The continuity equation of the dust of a certain grain 
size including radial drift and turbulent mixing reads 



Ek + -dr (rFk) = 
r 

The dust mass flux Fk is given by 

Fk = SkWd°*st.k - Dd,k^gdr ( — ^ 



(22) 



(23) 



The flrst term on the right side is the mass flux for the 
radial drift of individual particles discussed in the last 
section. The second term is the mass flux due to turbulent 
diffusion (radial mixing). 

2.6. Dust coagulation 



Two 
disk 



particles of mass rrii and rrij in a protostellar 
tend to hav e diff erent particle velocities Vi and 
20001 ) ■ For example, micrometer-sized 



Vj (jBeckwith et al 



particles are carried away with the gas while larger boul- 
ders, like meter-sized bodies, which are decoupled from 
the gas, are hardly affected by any gas motion. Small par- 
ticles have high relative velocities due to Brownian mo- 
tion even for equal-sized pairs. Larger particles have not. 
These relative velocities Awy lead to occasional collisions. 
The number of collisions per second between two particle 
species with number densities A^i and Nj can be calculated 

to be 

collisions . ,^ ,^ ,„ ,, 
^ = A^yayiViiVj . (24) 

We make the simplification that we take for Awy the av- 
erage relative velocity between the two dust particles. In 
principal, the relative velocities have stochastic variations, 
but we ignore them here. We assume that collisions lead to 
coagulation with a certain sticking probability Pc • In gen- 
eral this probability depends on various particle parame- 
ters like the size of the part icle, solid particle den sity, the 
'fluffiness' of the particles (JBlum &: Wurmll2000l ). It also 
depends on the relative particle velocity Avij. Small par- 
ticles tend to stick to each othe r up to high relative veloc- 
ities ( Dominik fc Tielenslll997l ) while larger bodies show 



the tendenc y to fragment even for small relative velocities 
( BenauOOOl) . This sticking probability will be discussed in 
more detail in section (|2.8p . 

With the collision rate Eq. ([24|) we can calculate the 
number of dust particles per second with mass mi which 
coagulate with any dust particle of any mass. 



N, 



Loss 



J2^Vii<^iiPcNiNj 



(25) 



which corresponds to the loss term for the number density 
of particles with mass ttij. The factor AvijaijPc is often 
called the coagulation kernel. The gain term for particles 
with mass TOj due to the coagulation of smaller particles 
with mass ?7ik and ttij reads 



N 



Gain 



7ni—m]^-\-mj 



AwkjCTkjPc^k^j 



(26) 



Now we can introduce the full co agulation equation which 
is given by (|Smoluchowski|[l91a) 



Ar_ jy-Gain tw^Loss 



(27) 



The continuous formulation of t he coagulation equation, 
as discussed by ISafronovl (|l969l ). corresponds to a non- 
linear integro-differential equation. However, this equation 
is rather difficult to solve both in its discrete or continu- 
ous version. This is only within the realms of possibility 
for very simple (and unfortunately unphysical) kernels. 
Therefore, we will solve the coagulation equation numer- 
ically. The algorithm we will make use of is described in 
detail in Appendix El 

2.6.1. Relative dust particle velocities 

We will consider four different sources for relative parti- 
cle velocities which lead to coagulation: Brownian motion, 
differential settling, turbulence and radial drift. 

First, let us focus on Brownian motion. Two particles 
of mass mi and m2 in a region of the disk with tempera- 
ture T have an average statistical relative velocity due to 
Brownian motion given by 



Awb = 



'SkT{mi + TO2) 

TTTO 17712 



(28) 



This expression shows that relative thermal velocities are 
higher for smaller dust particles than for larger dust parti- 
cles. Hence, the growth process due to Brownian motion is 
more effective for small particles than for large particles. 
For example, if we assume a temperature of 200 K, a solid 
particle density of 1 g/cm^ and micrometer-sized particles 
then the relative particle velocity due to Brownian motion 
is 0.2 cm/s. Particles of centimeter in size lead to a relative 
velocity of 10^^ cm/s. This particular example shows that 
there is practically no coagulation due to Brownian motion 
for particles much larger than micrometer size. In general, 
growth by Brownian motion leads to fractal str uctures and 
'fluffy' aggregates (cluster-cluster aggregates) (jOssenkopl 
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1993t iKempf et al.lll9 99'). However, we will ignore these 
intrinsic properties of the dust particles in the course of 
this paper an d assume a constant sol id material density . 
See, however. [Schmitt et all (|l997l ) or'Orm el et"al] (|2007t ) 
for dust particle coagulation models including porosity at 
a fixed radius in the disk. 

Differential settling is the second process that leads 
to relative velocities. If we assume that the solid par- 
ticles are smaller than the mean free path of the gas 
then the equilibri um settling velocity is given by zStf^k 
(JDullemond fc Dominik 2004 ). In this expression St is the 
Stokes number introduced in section 12.41 However, for 
Stokes numbers larger than unity, the equilibrium settling 
velocity model loses validity. Very large bodies (St -^ oo) 
above or below the midplane follow an orbit that is tilted 
with respect to the midplane. The settling velocity to- 
wards the midplane can not exceed the vertically projected 
Kepler velocity zVy^/r corresponding to this inclined orbit. 
For this reason we restrict the settling velocity to be the 
projected Kepler velocity at most. Considering this, we 
adopt the following settling velocity in our model. 



vs 



zStilk 
1 + St 



(29) 



The relative settling velocity between two particles of mass 
TOi and TOj at a height z above the midplane then reads 



Avg = zQi- 



St; 



Sti 



1 + Sti 1 + Stj 



(30) 



The third source for relative velocities of particles in 
the disk is the radial drift which was discussed in detail 
in section [2751 The relative velocity in this case is simply 
the difference in the drift velocities 



^«D = l^dust.i 



^dustj I 



(31) 



The fourth relative velocity between the particles 
is due to turbulence in the disk. Relative particle ve- 
locities pr oduced by t urbule nce w ere calculated i iumer - 



R 



ically by 'Volk et alj (|l980[ ) and iMizuno et aD (|l988l ). 
Weidensc hilling (1984!) fitte d these results w i th an alytical 
formulas. Current work bv lOrmel fc Cuzzil ( 20071 ) shows 
that these expressions underestimate the turbulent rel- 
ative velocities for particles with large Stokes numbers. 
In this paper, we wil l use the expressions calculated by 
Ormelfc Cu"izil (|2007l) . 

To give an impression of relative dust particle velocities 
in a protostellar disk. Fig. [3] shows a contour plot of this 
quantity at 1 AU including Brownian motion, differential 
settling, relative turbulent velocities and relative particle 
drift velocities. The same calculation at 10 AU is shown 
in Fig. m 

2.7. Particle fragmentation 

Collisions between particle aggregates do not necessarily 
lead to particle growth. For sufficiently high relative colli- 
sion velocities, the aggregates may fragment into smaller 
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Fig. 3. Relative velocities of dust particles at 1 AU in the 
disk as discussed in Sec. 12.6.11 This calculation includes 
Brownian motion, differential settling and relative turbu- 
lent velocities. In this calculation we adopted a turbulent 
a value of 10"'^. 
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Fig. 4. As Fig. [31 but now at 10 AU in the disk. 



bodies. The critical threshold velocity for this destructive 
process generally depends on the mass of the colliding 
particles. More precisely, fragmentation tends to play a 
non-negligible role if the kinetic collisional energy of the 
particles is of the order of the internal binding energy 
of the particles ( Borkowski &: Dweklll995h . Fragmentation 
velocities of aggregates are usually of the order of a few 
cm/s up to several 10 m/s. While smaller particles tend 
to st ick to each other up to h igh relative particle veloci- 
ties ( Dominik &: Tielenslll997l ) larger bodies show th e ten- 



dency to fragment even for small relative velocities ( Bena 
I2OOOI ) . For simplicity, we will assume a fixed threshold ve- 
locity for particle destruction Vf which does not depend 
on the mass of the particles. However, we will investigate 
how the results of the simulations change if V{ is varied 
over a wide parameter range. The dependency of Vf on 
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the particle mass will be investigated in the near future 
including laboratory results of dust particle collisions. 

The result of destructive collisions between solid par- 
ticles, i.e. the exact particle distribution after fragmen- 
tation, is still a matter of debate. Usually this particle 
distribution is described by a power-law, 

n{rn) dm ex m~^dm . (32) 

In this expression n{m)dm is the number of particles 
per unit volume within the mass range [m, m + dm] . 
The last decades involved vario us attempts to deter mine 
the fragmentation param eter ^. iMathis et al.l ( 19771) and 
also lDraine fc Led ( 19841 ) showed that the extinction and 
scattering of starlight by interstellar dust can be repro- 
duced by a power-law dependency following ^ = 1.83. 
Experimental studies found values for S, ranging between 
1.3 (low- velocity impacts) and 2 (catastro phic impacts) 
( Davis fc Rvanlll990l:lBrum fc Muenchlll993[ ). Steady state 
solutions between coag ulation and fragm entation lead to 
£ = 1.83 as s hown by iDohnanvil (19691). More recently, 
iTanaka et al.l (19961) argued that the very general result 
£ = 1.83 is a direct implication of the self-similarity of the 
particle size distribution. In this paper we will assume the 
£- value 1.83 if not otherwise noted. 

The process of fragmentation between particles which 
have the same mass is different from the fragmentation of 
particles whose mass differ by orders of magnitude. Two 
bodies of equal mass may destroy each other. Small dust 
grains, however, are not able to destruct a meter-sized 
body. But they can excavate a small crater in the larger 
target. This process is usually called 'cratering'. We will 
assume that cratering sets in if the mass of the colliding 
bodies differs by more than one order of magnitude. In this 
case, the smaller dust particle ms excavates a crater which 
contains a factor x times its own mass, i.e. mcratcr = X"^s- 
The parameter x is set to unity if not otherwise noted. The 
mass of the smaller body and the crater ejecta are then 
redistributed according to Eq. ([5^ . On the other hand, 
if the mass of the colliding particles differs by less than 
an order of magnitude, i.e. in the non-cratering case, then 
the total mass is redistributed following Eq. ([32|) . 

To illustrate the results of fragmentation Fig. [5] shows 
the outcome of a destructive collision as modeled in the 
paper at hand. The solid line shows the outcome of frag- 
mentation in the case of cratering. The dotted line corre- 
sponds to the fragmentation results of two particles with 
the same mass. 

2.8. Coagulation and fragmentation probabilities 

If the collision velocity of two boulders is sufficiently large 
then the particles tend to fragment into smaller bodies 
instead of coagulating to larger aggregates. We will assume 
that the probability for fragmentation pf only depends on 
the relative particle velocity Av and adopt the following 
expression for this probability, 

Vf 




pi{Av) 



9 {v[ — v) + Q {v — V[) 



(33) 
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Fig. 5. The assumed fragmentation results of coUisional 
destruction as discussed in Sec. 12.71 The solid line shows 
an example for cratering. The larger body has a mass of 
~ 10^ g. The smaller dust grain, which is destroyed in 
this process, had a mass of ~ 10~^ g. The dotted line cor- 
responds to the collision of two particles with the same 
mass. In this calculation the fragmentation parameter £ 
is assumed to be 1.83 which means that most of the frag- 
mentation results are at the large end of the particle size 
distribution. 



The two Heaviside step functions Q ensure that the par- 
ticles fragment with 100% probability if the relative par- 
ticle velocity Av is larger than the critical fragmentation 
velocity V{. For Av < Vf we assume that there is always 
a possibility for fragmentation given by (v/vi)'^. We will 
investigate the influence of the critical fragmentation ve- 
locity Vf and the index ip. The value of ip is set to unity 
if not otherwise noted. The probability for coagulation pc 
is given by Pc = 1 — Pf • The last expression implies that 
the particles either coagulate or they fragment. We do not 
allow the particles to collide and not to undergo either the 
process of coagulation or fragmentation. However, just for 
the moment let us assume that Pi + Pc < 1 • If this last 
expression holds then the time scales for coagulation and 
fragmentation increase. However, this issue might be con- 
sidered in a forthcoming paper. 

3. Simulation results 

In the following we will present the results of our numerical 
simulations. These simulations include various effects, for 
example, different particle growth mechanisms (Brownian 
motion, differential settling, etc.), the radial drift of the 
dust and particle fragmentation. To illustrate the influence 
of these effects on the particle growth, we will proceed in 
certain steps. In every step more and more effects will be 
included. In the first step we will consider the growth of 
the dust particles at various radii in the disk. In this part 
we do not allow the particles to move radially. In the sec- 
ond step, however, we will also include the radial motion 
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of the dust which was discussed in detail in section [^31 In 
the last step we will additionally consider fragmentation. 

3.1. Step 1 - Coagulation only 

What are the growth time scales of the solid particles 
at different radii in the disk? To answer this question, 
we will not allow any radial motion of the particles. We 
glue the dust to a certain radial position even though 
the radial drift of the dust is potentially very high. We 
also do not allow particle fragmentation. The coagula- 
tion of solid particles at a fi xed radius in t h e disk was 
for example also treat e d by ISchmitt et al.l (119971) an d 



(IMlD, 



Dullemond fc DominikI ( 20051 1 . iNakaeawa et ah. .^-_ 
Tanaka et all (|2005l ) and recentlv ICieslal (|2007t l. We as- 
sume that the mass of the disk is 1% of the central mass, 
an initial dust-to-gas ratio of eo ~ 10~^ and a solid mate- 
rial densitjoof ps = 1.6 g/cm'^. The turbulent a parameter 
is 10~^ and the turbulent q parameter is set to be 1/2. At 
the beginning of each simulation the dust is equally dis- 
tributed between a dust particle size of 0.5 /im and 0.8 
/xm. 

Let us first focus on the particle growth due to 
Brownian motion at different radii in the disk. The re- 
sult of this simulation, i.e. the particle size distribution 
after 1 Myrs, is shown in the upper panel of Fig. [6] 
According to these results dust particles grow from sub- 
micrometer to ~ 30 /im in radius in 1 Myrs at 1 AU in the 
disk. At 10 AU the particle distribution has a maximum 
for a ^ 4 /im. At 100 AU most of the dust is roughly 
a micrometer in size. We conclude that particle growth 
due to Brownian motion is not very effective, which is a 
well known result (IOssenkopJll993t ISchmitt et al.lll997 : 
Dullemond fc Domimkll2005r ). However, Brownian motion 
is an important effect for the following reason. We calcu- 
late the relative velocities due to Browian motion, differ- 
ential settling and turbulence for a = 0.6 /im equal-sized 
particles at 1 AU in the disk. While the relative parti- 
cle velocity due to Brownian motion is 0.4 cm/s the rela- 
tive turbulent velocity is in the order of lO^'* cm/s. The 
relative velocity due to differential settling is practically 
zero. Dust particle growth due to differential settling or 
turbulence gets of importance only for larger particles. 
Therefore, Brownian motion is a trigger mechanism for 
the entire coagulation process which was noted before by 
WeidenschiliinS (|l984l ). 

Now, we will additionally include coagulation due to 
differential settling into our model. The result of this sim- 
ulation is shown in the second panel of Fig. [6l This plot 
shows that particles have grown to more than 10"* cm in 
radius at 1 AU in the disk after 1 Myrs. This particle size 
is more than 6 orders of magnitude larger than the grain 
size after 1 Myrs caused by Brownian motion. Most par- 
ticles at 10 AU and 100 AU have grown to sizes of about 
1 m and 100 /im, respectively. We conclude that differen- 
tial settling is an effective growth mechanism which can 
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Fig. 6. These plots show the particle size distribution at 
different radii in the disk after 1 Myrs of disk evolution 
as discussed in Sec. 13.11 The left and the right plot al- 
ways belong together. From top to bottom more and more 
growth mechanisms are included in the simulations. The 
upper panel shows coagulation only due to Brownian mo- 
tion (BM). The second panel shows BM and differential 
settling (DS). Finally, BM + DS and turbulent coagula- 
tion (TC) are shown in the lowest panel. The left plots 
show the surface density of the dust at 1, 10 and 100 AU 
as a function of particle radius. On the right side the cor- 
responding contour plots of the dust surface density are 
shown as a function of the radial location in the disk and 
the particle radius. In these simulations the radial drift 
as well as the fragmentation of the dust particles were 
neglected. 



create large boulders in the inner parts of the disk. Note 
that in our model the vertical mixing continuously allows 
the grains to go back up again and grow again by dif- 
ferential settling. Therefore, the maximal size formula of 
Safronovl (J1969I ) does not apply here. 

Apart from the fact that particles grow to much larger 
sizes if differential settling is included. Fig. [S] also shows 
that there is always a certain amount of small particles 
that remains in the disk and that does not coagulate for 
at least 1 Myrs. After this time, roughly 6% of the dust 
between 1 and 75 AU is still present in grains < 1 mm. 
The reason for this is the following. Not all of the dust par- 
ticles coagulate at the same time. While a certain fraction 
of the dust has already grown to larger sizes and formed a 
thinner dust layer according to Eq. p3|) , a certain fraction 
of small dust remains in the higher regions above the mid- 
plane. These small dust particles high above the midplane 
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are subject to a rather slow coagulation process. The dust 
densities above the midplane are low after most of the dust 
already settled closer to the midplane. This leads to long 
growth time scales according to Eq. p4|) . Larger particles 
close to the midplane can not sweep up the smaller par- 
ticles above the midplane since turbulence is not able to 
stir them up so far. For this reason small particles remain 
in the disk for a long time. 

We will now also include relative velocities of the par- 
ticles caused by random turbulent motions. The result of 
this simulation is shown in Fig. [6] in the lower panel. This 
plot indicates that the dominant grain size at 1 AU, i.e. the 
grain size corresponding to the surface density maximum, 
changes by a factor of ~ 10 if turbulent coagulation is in- 
cluded in the simulation. The dominant particle radius at 
1 AU is ~ 10^ cm. At 100 AU, random turbulent motions 
also speed up the coagulation process which leads to par- 
ticles of a few centimeters in radius after 1 Myrs of disk 
evolution. Without relative turbulent velocities included 
in the simulation, the particle radius was two orders of 
magnitude smaller. 

3.2. Step 2 - Coagulation and radial motion 

We will now include radial motion, both as transport and 
as extra source of relative velocities for coagulation. This 
significantly changes the results of the last section. We 
find that the radial drift of solid particles is so high that 
the dust drift into the evaporation zone long before larger 
particles in the disk can possibly form. This happens even 
though an additional source for coagulation is introduced 
which decreases the coagulation time scales. We will in- 
vestigate if particles can in some way " break through" the 
radial drift barrier. 



3.2.1. Time evolution of the disk 

Fig. [7] shows the time evolution of the model. This plot 
indicates that cm-dm-sized particles form in the inner 
regions of the disk (< 2 AU) within the first 10'^ yrs. 
Compared to the outer parts of the disk the formation 
of these particles appears rather quickly due to compara- 
tively high gas and dust densities and high temperatures. 
With increasing distance from the central star the forma- 
tion of larger particles gets more and more difficult. At 
10 AU in the disk, it is still possible to form mm-sized 
particles in 10''' yrs according to Fig. [71 However, in the 
outer parts of the disk (> 100 AU) the dominant parti- 
cle size of the dust never exceeds 0.1 mm at any time. 
The disregard of radial drift in the previous section led to 
particle sizes of more than a centimeter at 100 AU after 
1 Myrs, which is orders of magnitude larger. 

The neglect of radial drift, as discussed in section [XTl 
involved a permanent amount of small particles which was 
present throughout the disk for at least 1 Myrs. These 
small particles were located high above the midplane and 
were subject to a rather slow coagulation process due to 



relatively low dust densities. Fig. [7] indicates that there is 
a smaller remaining amount of small dust if radial motion 
is taken into account. This is due to the following reason. 
Even the small particles in the higher regions of the disk 
can have relative radial velocities of the order of some 
mm/s or even cm/s. These higher relative velocites lead 
to higher collision rates and, hence, to a depletion of the 
small dust grains. 

After 10^ yrs of disk evolution, the average particle 
size at a certain radius in the disk starts to decrease in 
time. To give an example, after 10^ yrs the dominant dust 
grain radius at 1 AU is ^ 1 cm. After 1 Myrs this value is 
about an order of magnitude lower. While particles drift 
inward from a certain radial position they are replaced 
by other particles from the outer parts of the disk. The 
coagulation time scales are larger in the outer parts of 
the disk which means that particles grow to smaller sizes 
in the same time. Therefore, the particles that reach a 
certain position are smaller than the particles that drift 
away and, hence, the average dust particle size decreases. 

In the simulation shown in Fig. [7] the Stokes number 
of the dominant particles never exceeds unity, i.e. never 
breaks the radial drift barrier, at any disk radius consid- 
ered at any time. This is indicated by the St = 1 line 
which is also shown in this plot. At '--^ 1 AU in the disk, 
the simulation shows that particles may grow to sizes that 
correspond to a Stokes number slighly smaller than unity. 
In the following we will investigate if the particles may 
break through the St = 1 barrier for certain disk parame- 
ters. 

3.2.2. Effect of disk mass 

We investigate the effect of disk mass on the particle 
growth. The result of this investigation can be seen in 
Fig. [HI This plot shows the dominant dust particle size 
for different disk masses after 10"' yrs of disk evolution as 
a function of disk radius. We find that the particle size 
increases by an order of magnitude if the disk mass is in- 
creased from 1% to 20% of the central mass. Larger disk 
masses lead to higher gas and dust densities and, hence, 
to higher collision rates according to Eq. (|24|) . Therefore, 
dust particles can grow to larger sizes over the same time 
interval. 

The Stokes number of the dominant particles is always 
smaller than unity. Of course, particles may grow to larger 
sizes which increases the Stokes number since St oc a. 
However, larger disk masses also lead to higher gas den- 
sities which again decreases the Stokes number because 
St oc 1/pg. Finally, both effects cancel out and the disk 
mass seems to plays a minor role in breaking the radial 
drift barrier. 



3.2.3. Effect of turbulence 

As in the last section, we calculate the dominant parti- 
cle size after 10"* years but now for different turbulent 
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Fig. 7. The particle size distribution at different radii in 
the disk at different times of disk evolution as discussed 
in Sec. 13.2.11 In this simulation all particle gowth mecha- 
nisms are included as well as the radial motion of the dust. 
The fragmentation of particles is neglected. The left and 
the right plots always belong together. The left column 
shows the surface density as a function of particle radius 
at 1, 10 and 100 AU. The right column shows the corre- 
sponding contour plots of the surface density as a function 
of disk radius and particle radius. The white lines in the 
contour plots denotes the particle radius for which the 
Stokes number is unity (i.e. largest radial drift and largest 
radial velocities). 



a-parameters instead of different disk masses. The initial 
dust-to-gas ratio in this simulation is 10~^, the disk mass 
is 10~^ Mi, and the result is shown in Fig. [9l 

One would intuitively think that in a certain time par- 
ticles can grow to larger sizes in highly turbulent disks 
than in low-turbulent disks. Fig. [5] shows, however, that 
the dominant particle size after 10^ yrs is only weakly de- 
pendent on the turbulence parameter a. If a changes by 
two orders of magnitude then the dominant particle size 
only changes by a factor of two. This can be understood 
by the following consideration. 

A high amount of turbulence in the disk leads to high 
relative turbulent pa,rticle velocities (|V51k et al.l Il980t 
Weidenschilling|ll984l : [Cuzzi et al.ll200l[ ). These high rela- 
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Fig. 8. The effect of disk mass on the particle growth as 
discussed in Sec. 13.2.21 Shown is the dominant dust par- 
ticles radius after 10^ yrs of disk evolution for different 
disk masses between 0.2 and 10 AU. The turbulent a pa- 
rameter is 10~^ and the initial dust-to-gas ratio is 10^^. 
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Fig. 9. Same plot as Fig. [8] but now showing the ef- 
fect of turbulence on the particle growth as discussed in 
Sec. 13.2.31 Shown is the dominant particles size after 10^ 
yrs of disk evolution for different turbulent a parameters 
between 0.2 and 10 AU. The disk mass is 10"^ M^ and 
the initial dust-to-gas ratio is 10^^. 



five velocities cause high collision rates, cf. Eq. (p4)) . which 
favour the process of coagulation. For this reason parti- 
cles should have grown to larger sizes in highly turbulent 
disks. On the other hand, a large amount of turbulence 
leads to thick particle layers. The dust is stirred up in the 
higher regions of the disk which causes the average dust 
densities to decrease. The collision rates in Eq. p4)) are 
proportional to the particle number densities. Lower dust 
particle densities lead to longer coagulation time scales. 
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The two determining factors for the growth time scales, 
the relative turbulent particle velocity and the dust den- 
sity, seem to cancel out if the amount of turbulence in the 
disk is varied. Hence, different a-parameters lead to the 
same particle size over the same time interval. 

3.2.4. Effect of the initial dust-to-gas ratio 

We now investigate the effect of the initial dust-to-gas ra- 
tio on the growth time scales and the particle size distri- 
bution. We consider a disk mass of f 0~^ M^, and a turbu- 
lence parameter a of f 0^'^. The result of this investigation 
is shown in Fig. [TUl 

This contour plot shows the surface density of the par- 
ticle distribution as a function of disk location and par- 
ticle radius for four different initial dust-to-gas ratios af- 
ter 10^ yrs of disk evolution. These results indicate that 
f 0"* — f 0^ cm sized boulders can form in the inner parts of 
the disk (< 3 AU) subject to the condition that the initial 
dust-to-gas ratio of the disk is higher than f %. This means 
that the dust particles may overcome the radial drift bar- 
rier if the dust-to-gas ratio is slighly higher than usually 
assumed. A contour plot of the surface density distribu- 
tion with eo = 0.03, i.e. in the case where the particles are 
able to break through the radial drift barrier for disk radii 
< 3 AU, as a function of time is shown in Fig. [TT] 

To understand this importance of the initial dust-to- 
gas ratio we co nsider the growth rate of the dust particles 
as given bv iKornet et al.l (|2001l ). 
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The dust mass density can be approximated by pd oc Yi^/h 
so that we obtain 



a — — vaSt Cs 



(36) 



With the height of the dust layer Eq. (fT5|) . the last ex- 
pression can be written as (for St < 1) 

1 



-eo^f Stsli^ 



(37) 



If we also take into account the definition of the Stokes 
number in Eq. (fT3|) . then most quantities cancel each other 
out, particularly the gas surface density Eg, leading to 



a = aeoilk 



with the solution 



aoe 



eo^k^ 



(38) 



(39) 



This expression shows that only the initial dust-to-gas ra- 
tio eo and the Kepler frequency Slk determine the turbu- 
lent growth time scales as long as St < I. According to 
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Fig. 10. These plots shows the effect of the initial dust-to- 
gas ratio on the particle growth as discussed in Sec. 13.2.41 
The right side shows contour plots of the surface density as 
a function of disk location and particle radius for 4 differ- 
ent initial dust-to-gas ratios after 10^ yrs of disk evolution. 
The corresponding left plots show the surface density as 
a function of particle radius for 3 different locations in 
the disk (0.3 AU - sohd, 1 AU - dotted, 3 AU - dashed) 
after the same time. The disk mass is 10~^ M^ and the 
turbulent a parameter is I0~'^. For initial dust-to-gas ra- 
tios which are slighly higher than 1% the particles break 
through the radial drift barrier. 



Eq. pg]) . the time scales are not linear dependent on the 
initial dust-to-gas ratio. An increase of eo leads to an ex- 
ponential decrease of the growth time scales. This strong 
dependency unveils the crucial importance of this initial 
parameter. Eq. ([55)1 also shows that turbulent coagulation 
occurs faster in the inner parts of the disk than in the outer 
parts since Vt-k oc r~^-^ . For this reason, the particles first 
break through the radial drift barrier in the inner parts of 
the disk (cf. Fig.[lO|). 

In Section 13.2.21 and 13.2.31 we have seen that the dom- 
inant particle size only shows a weak dependency on the 
disk mass and the amount of turbulence in the disk. This 
can also be explained by Eq. ((37)) . The turbulent growth 
rate of the dust is neither dependent on the disk mass nor 
on the turbulent a parameter. Moreover, this expression 
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Fig. 11. This plot shows the results of a simulation in 
which the particles can break through the radial drift bar- 
rier as discussed in Sec. 13.2.41 Particle fragmentation is 
neglected in this simulation. Shown is the surface density 
distribution for the first 10"* yrs of disk evolution for an 
initial dust-to-gas ratio of 0.03 as a function of disk radius 
and particle radius. The disk mass is 10~^ Mi, and the 
turbulent a parameter is 10"'^. The right side is a contour 
plot of the surface density. The left side shows the abso- 
lute values of the surface density for 3 different disk radii 
(solid - 0.3 AU, dotted - 1 AU, dashed - 3 AU). 



also indicates that the disk temperature and intrinsic par- 
ticle properties Uke sohd density are rather unimportant as 
long as the Stokes number of the particles is smaller than 
unity and turbulence is the leading process that triggers 
coagulation. 

()2007l) 



However, Ormel et al 



have shown that the 
porosity of dust particles actually matters in the early 
phases of disk evolution. This discrepancy is due to 
the fact tha t the Eg. ([39]) only holds if St > a while 
Ormel et al.l (J2007l ) considered particles with St < a. 
Moreover, Brownian motion is the main source for rela- 
tive particle velocities for small dust grains in the early 
disk evolution while the derivation of Eq. ((39)) assumes 
that turbulence is the major source for relative particle 
velocities. 



3.2.5. The radial drift barrier 

Now we will estimate in which regions of the disk and 
under which conditions the solid particles can theoretically 
overcome the radial drift barrier. 

In section [3.2.41 we have seen that particle coagulation 
due to turbulence in the disk can be described by d = 
afi^eo . We define a particle growth time scale Tg by 



a 

1- 
a 



1 



^kEo 



(40) 



The parameter 7 measures how much the solid particle 
has to grow to cross the particle size region of fast radial 
drift, i.e. to overcome the radial drift barrier. We assume 
this parameter to have a certain value determined by the 
disk model and to be a constant throughout the disk. The 
largest radial drift velocity in the disk is approximately 
given by c^/V\^. We define a radial drift time scale Td by 



Td 



The ratio between these two time scales is given by 

2 



Td 



1_ 



(41) 



(42) 



In the last step we made use of Eq. (|15p . Now, the particles 
may overcome the radial drift barrier if the ratio Tg/rd is 
smaller than unity, i.e. if the growth time scales are smaller 
than the radial drift time scales. The parameter 7 is still 
indefinite. 

To specify this parameter we consider Fig. [10] in 
Section [3. 2. 41 These simulation results show for which ini- 
tial dust-to-gas ratio eo the particles break through the 
meter size barrier at a certain radius in the disk. We chose 
the parameter 7 in a way that the condition Td > Tg is in 
agreement with the results shown in this figure. This leads 
to 7 « 12. With this value, the particles should overcome 
the radial drift barrier if the inequality 



eo>12 



H 



(43) 



holds. 

The particles, which break through the radial drift bar- 
rier in Fig. I10[ have already drifted inwards. For this rea- 
son, the critical value given by Eq. [43] indicates the initial 
dust-to-gas ratio for which the particles most likely break 
through the radial drift barrier. The sufficient eo-value to 
overcome the radial drift barrier is presumably even lower 
than this value. 



3.2.6. Dust mass loss in the disk 

When particles drift into the evaporation zone, they are 
lost for the process of planetesimal formation. Hence, the 
question of how much solid material is actually lost due to 
its drift into the inner regions is of essential importance. 
We calculate the mass which is present in small (St < 1) 
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Fig. 12. The mass of the dust disk for small (St < 1) and 
large (St > 1) particles between 0.5 AU and 150 AU as 
a function of time for different initial dust-to-gas ratios 
as discussed in Sec. 13.2.61 In this simulation only the par- 
ticle coagulation and the radial motion of the dust were 
considered. Particle fragmentation was neglected. 



and large (St > 1) particles between 0.5 AU and 150 AU 
as a function of time for different initial dust-to-gas ratios. 
The result of this calculation can be seen in Fig. [121 

This plot shows, that the mass of the dust disk does not 
change significantly within the first 10^ yrs for every initial 
dust-to-gas ratio considered. Since the power law index of 
the surface density is —0.8, most of the solid particles are 
in the outer regions of the disk. Fast radial drift in the 
inner disk regions which takes place in ~ 10'^ yrs does not 
change the total dust mass in the disk. 

After a few 10"* yrs, the mass present in small grains 
starts to decrease. After 1 Myrs of disk evolution, this 
mass is less than 1% of the initial dust mass. The amount 
in small grains, i.e. St < 1 particles, is dependent on the 
initial dust-to-gas ratio. For ep = 0.01 roughly 0.4% of 
the initial particle mass is present in small grains. For 
eo = 0.03 this mass is a factor of 4 lower. Hence, higher 
initial dust-to-gas ratios lead to lower dust masses in small 
grains after 1 Myrs. In the last section we showed that 
particles grow faster with increasing dust-to-gas ratio. 
Therefore, particles can grow to larger sizes while moving 
radially inwards. However, larger sizes also lead to higher 
radial drift velocities (Eq. [T7)) . For this reason, the mass 
which is present in small particles in the disk decreases 
faster for increasing dust-to-gas ratios. We will find the 
same behaviour for the small particles in step 3 where 
fragmentation is also taken into account. 

For an initial dust-to-gas ratio of 1%, the mass of the 
entire dust disk, i.e. the mass in small and large parti- 
cles, after 1 Myrs between 0.5 AU and 150 AU is 0.4% 
of the initial dust mass. Most of the dust has drifted into 
the evaporation zone. For higher initial dust-to-gas ratios, 
i.e. higher than 0.015, the particles in the inner regions of 



the disk can break through the radial drift barrier. These 
larger boulders around 1 AU then sweep up smaller par- 
ticles which drift inwards from larger radii (cf. Fig. [12] be- 
tween 10^ yrs and ^ 5 x 10^ yrs). After ^ 5 x 10^ yrs most 
of the dust mass is present in large boulders. While for 
eo = 0.015 roughly 20% of the initial dust mass is present 
in St > 1 particles after 1 Myrs, the remaining mass in 
large boulders is a factor of ~ 4 higher for eo = 0.03. 
Note that the mass of the remnant dust disk after 1 Myrs 
changes by a factor of ~ 200 by changing the initial dust- 
to-gas ratio from 1% to 3%. We conclude that the initial 
dust-to-gas ratio is a crucial parameter which has an im- 
portant influence on how much solid material remains in 
the disk after 1 Myrs. However, the mass present in small 
grains is always less than 0.4% of the initial dust mass 
after 1 Myrs no matter the value of e. 

3.3. Step 3 - Coagulation, radial motion and 
fragmentation 

We now also include particle fragmentation in our simula- 
tion. We investigate how this destructive effect influences 
the particle growth in the disk and how various disk pa- 
rameters, like the initial dust-to-gas ratio, the turbulence 
parameter a or the fragmentation velocity Vf, influence 
the coagulation/fragmentation process. 

3.3.1. Time evolution 

The evolution of the disk in the first 1 Myrs is shown 
in Fig. [T21 In this calculation, the fragmentation velocity 
is V[ — 10'^ cm/s and the fragmentation parameter ^ is 
1.83. We adopt a disk mass of 10"^ M^, a turbulent a- 
value of 10~^ and an initial dust-to-gas ratio of 10~^. The 
cratering-parameter x is 0.5 and V' = 2. 

After 10'^ yrs of disk evolution, most of the particles in 
the disk < 3 AU have grown to sizes of some millimeters. 
However, if fragmentation is neglected (cf. Sec. 13. 2p the 
dominant particle size at 1 AU in the disk after lO'^ yrs 
is an order of magnitude larger. This significant difference 
is due to the fragmentation of particles. When the par- 
ticles reach millimeter size then destructive effects pre- 
vent the particles from growing to larger sizes (cf. Fig. [3] 
with 10 m/s). Even after 10^ yrs, the dominant particle 
size in the disk < 10 AU is still of the order of a mil- 
limeter. Hence, this particle size corresponds to the frag- 
mentation barrier for this specific set of disk parameters. 
Even for long periods of time the particles are not able 
to overcome this barrier. Once the particles have reached 
the fragmentation barrier the particle distribution is char- 
acterised by an equilibrium between particle coagulation 
and particle fragmentation due to destructive collisions. 
In other words, the amount of particles of a certain mass, 
which are created by dust particle coagulation, equals the 
amount of particles, which are destroyed by high veloc- 
ity collisions. This steady state will be discussed in more 
detail later in this Section. 
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Fig. 13. As Fig. [71 but now also the fragmentation of 
particles is included in the simulations as discusssed in 
Sec. 13.3.11 The left column shows the surface density as a 
function of particle radius at 1, 10 and 100 AU. The right 
column shows the corresponding contour plots of the sur- 
face density as a function of disk radius and particle size. 



Fig. 1131 indicates that the maximum dominant particle 
size Omax and the Stokes number St have the same radial 
behaviour. This is due to the fact that relative particle ve- 
locities in our model (except Brownian motion) scale with 
this dimensionless number. For this reason, the dominant 
particle size follows Omax ex r~°'* which we obtain directly 
from the definition given by Eq. p3p . 

Due to destructive collisions a large amount of dust is 
present in small grains as can be clearly seen in Fig. [131 
We calculate the amount of dust which is present in grains 
larger (smaller) than 10~^ cm after 10^ yrs of disk evo- 
lution. While 18% of the dust mass is present in grains 
larger than 10~^ cm, yet 82% of the mass is present in 
smaller grains. This large population of sub- mm grains 
should have a strong effect on the spectrum of the proto- 
stellar disk. However, we will not investigate the influence 
of the fragmentation parameters, i.e. V{ and ^, on the disk 
spectrum which goes beyond the scope of this paper. This 
will be investigated in the near future. 
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Fig. 14. The influence of the turbulence parameter a on 
the dominant particle size after 10^ yrs of disk evolution 
for different disk radii between 1 and 20 AU as discussed in 
Sec. 13.3.21 The disk mass is 10~^ M^,, the fragmentation 
velocity is 10^ cm/s and the initial dust-to-gas ratio is 
10~^. This graph also shows the particle size for which 
the Stokes number ist unity. The x parameter is set to 0.5 
and -tjj — 1. 



3.3.2. Effect of turbulence 

Different turbulent a-values should lead to different max- 
imum particle sizes due to destructive collisions. To in- 
vestigate the influence of turbulence on the fragmentation 
barrier, we calculate the dominant particle size for differ- 
ent a- values after 10* yrs of disk evolution. In this simula- 
tion the disk mass is 10"'^ M*, the fragmentation velocity 
is 10"^ cm/s, the initial dust-to-gas ratio is 10~^ and the 
results of the calculation are shown in Fig. [HI 

According to this plot, the dominant particle size is 
fairly dependent on a in moderately turbulent disks. If a 
is changed from 10~^ to 10"'* then the dominant particle 
size Odom changes by a factor of '^ 5. We find that less 
turbulence shifts the fragmentation barrier towards larger 
particle sizes. Hence, in less turbulent disks particles can 
grow to larger sizes than in highly turbulent disks. 

However, this statement does not hold for extremely 
low turbulent disks. In these disks, turbulence is not the 
main source for relative velocities and, hence, the frag- 
mentation barrier should not be dependent on a. If a is 
smaller than ~ (cs/21^)^ (cf. Eqs. [HI and [321) which is 
~ 10""* at 1 AU then relative particle velocities due to ra- 
dial motion exceed relative dust particle motions induced 
by turbulence. To illustrate this independency we calcu- 
late the dominant particle size after 10'* yrs for a disk 
with a very low a-value of 10~^°. The result of this cal- 
culation is also shown in Fig. 1141 In this nearly laminar 
disk, destructive coUsions due to relative drift velocities up 
to 50 m/s prevent particle growth to sizes of more than 
-- 2 mm at 1 AU. 



16 



Brauer, Dullemond and Henning: Coagulation and fragmentation of grains 



100 




1 10 

Particle radius [cm] 



100 



-. 10 



100 









' ' 


■ 


^ 


0.005 


- 


r 


\ - 


E 






- ' ~y 


\ 


o 






- ' ' ' / * 


\ 


Ui 


0.000 




_^^ '> 


\ ■ 






w 








>> 


























1 




d 










(1) 










o 












-0.005 


- 


1 


' 


<> 










o 








1 






















-i 










{/) 








/ 




-0.010 


- 







0.005 



0.000 



-0.005 ^ 



-0.010 



0.01 



0.10 1.00 

Particle radius [cml 



10.00 



Fig. 15. The relative particle velocities at 1 AU as a func- 
tion of particle radius as discussed in Sec. 13. 3.^ The turbu- 
lence parameter a in this calculation is 10^^*^. This means 
that relative radial motion is the main source for rela- 
tive particle velocities. A critical fragmentation velocity 
of 10 m/s results in a very narrow band in which particle 
coagulation is still possible. If the particle size dispersion 
is larger than the extent of this bottleneck then particle 
fragmentation starts to play a non-negligible role. 



Fig. 16. The dust particle distribution (solid) and its time 
derivative (dashed) as a function of particle size after 
700 yrs of evolution at 1 AU in the disk as discussed in 
Sec. 13.3.21 The distribution is located around oq = 3 cm 
and it has a size dispersion of ~ 1.5 cm. In these 700 yrs, 
particle fragmentation was neglected. In the calculation of 
the source terms, which are shown in this figure, fragmen- 
tation is included. These source terms show that destruc- 
tive collsions would rapidly shift the dominant particle 
size to smaller values. 



Relative radial drift velocities are always due to par- 
ticle size differences. Monodisperse distributions do not 
show relative radial motion. The simulation result for ex- 
tremely low turbulent disks (a — 10~^°) raises the ques- 
tion how the particle size dispersion of the dust distribu- 
tion can produce such high relative velocities to inhibit 
particle growth to larger sizefl 

We try to answer this question by considering the rel- 
ative velocities of dust particles at 1 AU in the disk as 
a function of particle radius, c.f. Fig. [TS] In this calcu- 
lation, we adopt an a-value of 10"^*^ which means that 
relative radial motion is the main source for relative ve- 
locities. According to this figure, particle coagulation is 
only possible in a very narrow particle size interval, i.e. in 
the dark shaded regions of this plot. If the particle size dis- 
persion is larger than the extent of this 'bottleneck' then 
particle fragmentation starts to play a non-negligible role. 
With Eq. ([20|l for the radial velocities, we can estimate 
the importance of fragmentation for a specific particle size 
dispersion. We assume a particle size distribution which 
has a surface density maximum at oq = 3 cm. If the size 
dispersion is larger than 1 cm, then particles start to frag- 
ment with 100% probability. For a particle size dispersion 
of 0.5 cm and 0.1 cm, the fragmentation probability de- 
creases to 50% and 10%, respectively. Hence, only for par- 
ticle size dispersions of some millimeters, particles might 
have the chance to overcome the fragmentation barrier. 
For larger size dispersions, the fragmentation probability 



^ We define the particle size dispersion as the half-width of 
the size distribution 



is far too high to allow the distribution to pass the bot- 
tleneck shown in Fig. [151 

To investigate if the particle size dispersion is narrow 
enough to overcome the fragmentation barrier, we con- 
sider the following. We simulate 700 yrs of dust particle 
evolution neglecting fragmentation. The result of this sim- 
ulation, i.e. the particle distribution at 1 AU in the disk 
as a function of particle size, is shown in Fig. [TB] (solid 
line). The size dispersion of this particle distribution is 
^ 1.5 cm. Now, particle fragmentation tends to smear 
out the dust distribution and it increases the particle size 
dispersion. Hence, the distribution shown is Fig. [16] rep- 
resents a best case scenario; the distribution can not be- 
come narrower. What happens if we now switch on frag- 
mentation? Fig. [16] also shows the time derivative of the 
particle distribution (dashed line) if fragmentation is con- 
sidered. This curve indicates, that destructive collsions 
would rapidly shift the dominant particle size towards 
smaller values. The size dispersion is apparently too large 
for the particles to pass through the fragmentation bot- 
tleneck without undergoing substantial destructive coll- 
sions. If fragmentation is included in the simulations from 
the very beginning, then the size dispersion is even larger 
and, therefore, the chance of passing the narrow region of 
coagulation becomes even smaller. 

For a fragmentation velocity of 10 m/s, which we adopt 
in these simulations, the particles never overcome the frag- 
mentation barrier, regardless of the amount of turbulence 
in the disk since the radial drift always accounts for de- 
struction. We find that this statement also holds for larger 
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'(/'-values. We conclude that the amount of turbulence in 
the disk alone does not determine whether particles can 
break through the fragmentation barrier or not. Note, that 
the maximum radial drift velocity of particles is indepen- 
dent of radius, so that these statements hold everywhere 
in the disk. 

There are possible scenarios in which the radial parti- 
cle velocity, which is the main reason for particle frag- 
mentation in low turbulent disks, is lower than in the 
model discussed here. In these cases, particles might over- 
come the fragmentation barrier. One possibility are local 
gas pressure fluctuations. Since the radial drift velocity 
is proportional to the radial gas pressure gradient, local 
gas maxima can slow down and even prevent radial par- 
ticle motion. Therefore, we expect dust coagulation in- 
stead of dust fragmentation in these maxima. Also local 
dust particle enhance ments can slow down radial drift. 



Johansen et al 



( 20061 ) have shown that the radial drift 
velocity can be reduced by a factor of around 2. Further 
investigations of particle growth under these conditions, 
which make particle fragmentation less hkely, are needed. 
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Fig. 17. The dominant particle size as a function of 
disk location for 3 different fragmentation velocities af- 
ter 10^ yrs of disk evolution as discussed in Sec. 13.3.31 In 
this simulation, ^ — 2, x — 0, ^o — 0.03 and the turbulent 
a- value is 10~^. 



3.3.3. Effect of the fragmentation velocity 

For which critical fragmentation velocities can particles 
break through the fragmentation barrier? To answer this 
question let us consider a best case scenario. We adopt a 
low turbulent disk, i.e. a disk in which the relative radial 
velocities exceed the relative turbulent particle velocities, 
and we neglect the effect of cratering for the moment. We 
calculate the dominant dust particle size as a function of 
disk location for 3 different fragmentation velocities after 
10* yrs of disk evolution. The results of this calculation 
can be seen in Fig. 1171 In this simulation, the a-valuc is 
10~^, X = (no cratering) and tp ^ 2. 

For a fragmentation velocity of 5 m/s, particles can 
grow to millimeter size at ^ 1 AU in the disk before de- 
structive collisions prevent further particle growth. In the 
outer regions, i.e. at 10 AU, the dominant particle radius 
is a factor of 10 smaller. Even for a relatively high criti- 
cal velocity of 20 m/s the particles are not able to grow 
beyond a centimeter at 1 AU. 

For even higher fragmentation velocities, i.e. Vf ^ 
30 m/s, solid particles start to break through the frag- 
mentation barrier. Fig. [TH] shows the dust particle distri- 
bution for this critical velocity as a function of disk radius 
and particle radius for 4 different times of disk evolution. 
This plot indicates that particles have grown to meter size 
in the inner parts of the disk after 10* yrs. However, a 
fragmentation velocity of several 10 m/s for centimeter- 
or even meter-sized boulders is at least questionable. For 
lower (and probably also more realistic) critical velocities, 
i.e. velocities of 1 ... 10 m/s, we never find solid particles 
in our simulations which are able to overcome the frag- 
mentation barrier for any disk parameters considered. For 
a-parameters, which are higher than the adopted value of 
10~^ in the simulations of this paragraph, it is even more 



unlikely that solid particles may grow to larger sizes. This 
chance does not increase if destructive effects due to cra- 
tering are also taken into account. 

3.3.4. Disk dust mass 



As in Sec. 13.2.61 we calculate the solid material mass in 
the disk as a function of time, but now with the effect 
of particle fragmentation included in the simulations. The 
result of this calculation is shown in Fig. [191 

The dust mass does not change significantly within 
the first 10^ yrs for any eo considered. This is the same 
behaviour as in the case of no fragmentation. After 10^ yrs 
the mass starts to decrease rapidly. For an initial dust-to- 
gas ratio of 0.01 only 2% of the initial solid material mass 
between 1 and 150 AU remains after 1 Myrs. Higher initial 
dust-to-gas ratios lead to less solid material after 1 Myrs. 
For example, for e = 0.03 the mass is only 0.7% of the 
initial dust mass, which is a factor of ~ 3 lower. 

Let us compare the solid material mass after 1 Myrs for 
eo — 0.01 in the case of fragmentation/no fragmentation. 
We find that the remaining dust mass is a factor of 5 
higher if we allow the particles to destroy each other. This 
difference is due to destructive collisions which lead to 
large amounts of small particles in the disk (cf. Sec. 13. 3.1]) . 
These small dust grains have low radial drift velocities 
and, hence, long radial drift time scales. In other words, 
small particles stay much longer in the disk before they 
evaporate in the inner regions of the disk. For this reason, 
the solid material mass after 1 Myrs is higher in the case 
of fragmentation than in the case of no fragmentation. 

In the previous sections we found that if fragmentation 
is included in the simulations then the dust particles are 
not able to break through the meter size barrier. No larger 
particles in the inner parts of the disk can form which can 
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Fig. 18. These plots show how the particles break through 
the radial drift barrier and the fragmentation barrier as 
discussed in Sec. 13.3.31 Shown are contour plots of the 
surface density as a function of disk radius and particle 
radius at 4 different times of disk evolution. The fragmen- 
tation velocity is chosen to have the relatively high value 
of 30 m/s. In this simulation -0 = 2 and x = 0. The initial 
dust-to-gas ratio is 0.03. 



sweep up smaller dust particles drifting inward from the 
outer regions. For this reason, most of the solid material 
after 1 Myrs has drifted into the evaporation zone and is 
lost for the process of planetesimal formation. 

3.3.5. Effect of disk model 

In the introduction we mentioned that the disk model 
adopted in this paper differs significantly from the MMSN 
model. This leads to the question of how the results of this 
paper change if different disk models are considered. In 
this section, we repeat simulations of Sec. 13.31 with other 
disk model parameters, attempting to unveil the basic 
changes in the dust particle distribution. Table 1 shows 
the disk parameters for the simulations in this section. 
Model A and B are the MMSN model and the disk model 
in this paper, respectively. Model C is our model, but now 
with 10% disk mass instead of 1% compared to Mi,. This 
leads to gas densities which are comparable to those of 
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Fig. 19. The mass of the dust disk between 1 AU and 
150 AU as a function of time for 3 different initial dust- 
to-gas ratios as discussed in Sec. 13.3.41 In this simulation, 
particle growth particle fragmentation and radial motion 
are included. The initial disk mass of gas and solid mate- 
rial is 10-2 M*, a = lO"'*, X = 0.5, ■0 = 2 and vi = 10 m/s. 



Model 


Surface density 


Disk 


Temperature 




power law index 5 


mass 


power law index /3 


A 


1.5 


0.01 


0.50 


B 


0.8 


0.01 


0.50 


C 


0.8 


0.10 


0.50 


D 


1.5 


0.01 


0.62 


E 


0.8 


0.01 


0.62 


F 


0.8 


0.10 


0.62 



Table 1. Disk parameters for the simulations performed 
in Sec. 13.3.51 The quantity /3 denotes the temperature 
power law index T oc r"^. The Models A and B corre- 
spond to the MMSN model and the model adopted in this 
paper, respectively. Model C is as the model in this pa- 
per but now with 10% disk mass. The Model D to F are 
as A to C but with a slighly steeper radial temperature 
dependency. 



the MMSN model. The mass distribution, however, has 
a much flatter radial dependency. The models D to F are 
the same as A to C, but with a steeper ra dial temperature 
dependency. [Andrews fc Williams! ( 20071 ) observationally 
find radial temperature profiles with a median power law 
index of 0.62. This is slighly higher than the passively ir- 
radiated disk profile of 0.5 adopted in our model. 

Before we come to the results of the simulations, 
we will qualitatively discuss the difference between the 
MMSN model and the model in the paper at hand. The 
gas mass densities of our model are generally smaller than 
those of the MMSN model. This has the following main 
implications. First, solid particles are less coupled to the 
motions of the gas. The coupling between the gas and the 
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dust can be described by the Stokes number St, which 
is given by St = psafS. If the surface density of the 
gas S decreases, then St is shifted towards higher val- 
ues. Therefore, the particle growth barrier due to radial 
drift and particle fragmentation, which is usually referred 
to as the 'meter size barrier' and which corresponds to 
the particle radius implied by St = 1, is shifted towards 
lower particle radii. In the MMSN model, particles with a 
Stokes number of unity have radii of ~ 2 m at 1 AU in the 
disk. A surface density slope oi S = 0.8 implies a ~ 5 cm 
for St = 1 particles at 1 AU. While it seems challenging to 
grow particles larger than meter in size in the MMSN disk 
model, it is difficult to grow particles larger than centime- 
ter size in the disk model adopted in the paper at hand. 

Second, if the Stokes number is shifted towards higher 
values then all quantities depending on this number are 
influenced by this change as well. For example, for Stokes 
numbers smaller than unity the radial drift velocity of 
solid particles in t he disk is proportional to the Stokes 
number, Vr oc St (jWeidenschilling 1977a) . Now, if the 
Stokes number is modified due to a change of S then also 
the radial drift of the dust is significantly affected. The 
Stokes number also determines relative dust particle ve- 
locities in turbulent disks and, hence, dust particle growth 
time scales and the maximum dust particle size due to 
fragmentation. 

Fig. [20] shows the particle distribution after 1 Myrs of 
disk evolution for the Models A to F. In these simulations, 
particle growth, radial particle motion and destructive col- 
lisions are included. The initial dust-to-gas ratio is 10^^ 
and the a- value is 10""^. The -^-parameter is chosen to be 
2 and the cratering parameter x = 0-5. This figure shows 
that particles can grow to much larger sizes in model A 
than in model B in the inner parts of the disk. This is due 
to higher gas densities in the MMSN model which alter the 
Stokes number and shift the whole particle growth prob- 
lem towards larger particle radii. At 1 AU, the gas density 
in model A is a factor of ^ 15 higher than in model B. 
The dominant particle size before fragmentation inhibits 
further particle growth is 3 mm in model A and 0.2 mm 
in model B. This dominant particle size difference from 
one model to the other nicely mirrors the gas density dif- 
ference between the two models. Hence, we find that the 
dominant particle size is directly proportional to the gas 
density. 

Model C is the same as the model in our paper (B), 
but now with 10% disk mass instead of 1% compared 
to M^. Fig. [20] shows that the dominant particle radius 
due to destructive collisions is shifted by a factor of 10 
towards larger particle sizes. According to these results, 
particles can grow to a few millimeter in size in high 
mass disks before particle fragmentation prevents further 
growth. However, even in these very high mass disks, par- 
ticles can not overcome the fragmentation barrier. Since 
the whole coagulation/fragmentation process scales with 
gas density, higher disk masses do not provide a solu- 
tion for planetesimal formation. The entire particle growth 
problem is only shifted towards larger particle radii. 
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Fig. 20. The particle distribution in the disk after 1 Myrs 
for the disk models A to F as discussed in Sec. 13.3.51 The 
Models A and B correspond to the MMSN model and the 
model adopted in this paper, respectively. Model C is as 
the model in this paper but now with 10% disk mass. The 
Model D to F on the right side have a slighly steeper radial 
temperature dependency. 



The right column shows the results of the three simu- 
lations A-C if the radial temperature dependency follows 
T ex r^*^'^^ corresponding to the observational median. 
We do not find a significant difference in the maximum 
particle size between these two model sets. 

We also calculate the mass of the dust disk which is 
shown in Fig. [^TJ This plot shows that the remaining 
dust mass after 1 Myrs of disk evolution is smaller in the 
MMSN model than in the model adopted in this paper. 
This is due to the fact that the maximum radial drift 
velocity is proportional to the power law index S of the 
surface density profile (cf. Eg. [T8|). Since the parameter S 
is larger in the MMSN model than in our model, the max- 
imum radial drift speed is also larger. A higher drift speed 
leads to shorter drift time scales and, hence, reduces the 
remaining amount of dust after a certain time. 

In the disk models D-F, the temperature is generally 
smaller than in the models A-C. Therefore, the radial drift 
velocity is also smaller since «„ oc T. Hence, the disk dust 
mass in the model A-C after a certain time is generally 
smaller than in the models D-F. Finally, we find that less 
than 6% of the initial dust mass is left after 1 Myrs of disk 
evolution in any disk model considered. 
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Fig. 21. The mass of the dust disk as a function of time 
for the 6 different disk models A to F as discussed in 
Sec. 13.3.51 The models with 1% disk mass are normalised 
to unity. The disk models with 10% disk mass are nor- 
malised a factor of 10 larger. 
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Fig. 22. The effect of cratering on the equilibrium particle 
distribution as discussed in Sec. 13. 3.^ Shown is the surface 
density of solid particles at 1 AU in the disk after 10''' yrs of 
disk evolution for different cratering-parameters x- In this 
simulation, the fragmentation velocity is 20 m/s, e = 0.03, 
a = 10-5 and V = 2. 



3.3.6. Effect of cratering 

If a smaller particle of mass rris collides with a larger body 
at a sufRciently high velocity then the smaller particle does 
not only fragment due to this destructive collision but it 
can also excavate a certain amount of matter ttIc from the 
larger body, i.e. ttIc = X'^s- This effect is called cratering. 
In the following we investigate if this process has an effect 
on the equilibrium particle distribution between particle 
coagulation and particle fragmentation. 

Fig. [52] shows the equilibrium particle distribution at 
1 AU in the disk after 10^ yrs of disk evolution for differ- 
ent cratering-parameters x- The fragmentation velocity is 
20 m/s, e = 0.03, a = lO"""^ and ip = 2. 

This plot shows that the equilibrium distribution is 
hardly affected whether the effect of cratering is included 
in the simulations or not. Changing the x-value from (no 
cratering) to 1 (the projectile particle excavates a crater 
corresponding to its own mass) changes the surface density 
for 50 /im-sized dust grains by a factor of 1.3 at most. The 
maximum peak of the surface density is shifted from 9 mm 
to 6 mm by including cratering. We also investigated the 
effect of cratering for different fragmentation velocities. 
In any case we found that cratering does not significantly 
affect the particle distribution. This shows that the main 
destruction by fragmentation is due to collisions between 
particles of not large mass ratio. 

For laboratory experiments which investigate the colli- 
sion of particles for astronomical purposes it is interesting 
to know which particles collide with which particles in the 
equilibrium solution discussed in Sec. 13.3.11 For this rea- 
son, we calculate collision rates between dust particles in 
the disk in Appendix [E] and we show which collisions are 
most important for growth/fragmentation. 



4. Discussion and conclusions 

In this paper we study the evolution of the dust in 
a protoplanetary disk. W e expa nded on the work done 
by iDullemond fc DominikI (J2005l ) by adding an improved 
treatment of dust aggregate fragmentation and by in- 
cluding the self-consistent radial drift and mixing of dust 
aggregates. In addition to this, our integration method 
is faster by orders of magnitude than the original work 
by Dullemond & Dominik, by virtue of two new meth- 
ods. The first of these is the integration of the equations 
in vertical direction without deleting the vertical struc- 
ture information, which can be done because the verti- 
cal sedimentation-mixing equilibrium is assumed to be 
known at all times. Secondly, we integrate the coagulation- 
fragmentation in time using implicit integration, allowing 
us to overcome the extreme time scale stiffness that pre- 
vented Dullemond & Dominik from performing full scale 
models including fragmentation over time scales of mil- 
lions of years. 

Using these models we are able to study the problem 
of the radial drift and fragmentation barriers in full de- 
tail. Our models show that, consistent with current be- 
liefs, the combination of radial drift and fragmentation 
is a strong limitation to growth of aggregates in disks. 
Typically, aggregates cannot grow to sizes larger than 
millimeters throughout the disk if threshold fragmenta- 
tion velocities up to several m/s are considered, though 
the precise maximum size is disk model dependent. This 
statement holds regardless of the amount of turbulence in 
the disk. For highly turbulent disks, it is the turbulence- 
induced relative velocities that cause the damage to aggre- 
gates, while for nearly non-turbulent disks it is the differ- 
ential radial drift that limits the growth. Only if we set the 
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fragmentation threshold velocity to the unlikely value of 
30 m/s or larger do we find that the fragmentation barrier 
can be broken and particles grow to larger sizes. Whether 
this high fragmentation threshold is realistic remains to be 
verified by high-speed laboratory collision experiments. It 
has been tentatively shown that high-speed impacts of a 
small projectile on a large target may result in growth 
( Wurm et al.ll2005[ ). but further study in this direction is 
imperative. 

We have also investigated what happens when, for 
some reason, no fragmentation occurs. We then find that 
even in this reduced model, the particle radius never ex- 
ceeded several centimeters at any time at any radius in 
the disk because of radial particle motion. However, we 
demonstrated that this radial drift barrier problem is very 
sensitive to slight changes in the initial dust-to-gas ratio. If 
slighly higher initial dust-to-gas ratios than the canonical 
value of 1% are adopted in the simulations, then particles 
can grow to very large sizes in the inner parts of the disk. 

We do not include non-linear feedback of the dust back 
onto the gas in our rn odel. It has been recently shown by 
Johansen et al.l ( 20071 ) that such feedback can lead to the 
rapid formation of gravitationally bound clumps of dust 
which subsequently form Ceres-size bodies. The 'dust' par- 
ticles, however, must be large (Stokes number near unity) 
before this scenario can take place. We find that for low- 
turbulent disks Stokes numbers larger than 0.1 can be 
reached, but we need further investigation if the amount 
of dust present in these large grains is sufficient to trigger 
such a gravitational collapse in locally overdense regions 
in the midplane of the disk. 

In the full model, most of the solid material has drifted 
into the evaporation zone after 1 Myrs and the rem- 
nant disk contains less than 5% of the initial dust mass. 
However, there are reasons to believe that the strong ra- 
dial drift of particles in such disks may be reduced by 
non-linear hydr odynamic effects. Tentative results from 



Johansen et al.l |2006) find a reduction by a facto r of 3 
in MRI turbulence. Moreover, iBrauer et al.l (|2007l ) show 
that observations of millimeter grains in the outer regions 
of protoplanetary disks indicate that the standard ra- 
dial drift formulae are inconsistent with the observations. 
Hence, it may be important to investigate what happens 
to our model if radial drift velocities are reduced by some 
factor. This will be the topic of a fut ure paper. 

Th e d ynamic c alculatio ns o f iBarge fc Sommerial 
( 19951 ) and lKlahr fc H enning ( .1997[ ) also suggest that par- 
ticle trapping in pressure maxima could be a solution also 
to the relative velocity fragmentation barrier. We intend, 
in future work, to study how the dust evolves if such pres- 
sure maxima are present. 

By being able to model the dust evolution in a self- 
consistent way, our code may act as a basis upon wich 
further modeling of disks is done. For example, models 
of chemistry in disks are very dependent on the total 
surface area of dust grains per unit volume, in particu- 
lar when grain surface chemistry is taken i nto account 
( Aikawa fc Herbstlll999l : ISemenov et al.ll200d ). Models of 



MRI turbulence in disks depend strongly on the abun- 
dance of free electrons. This abundance is also very de- 
pend ent on the total surface area of dust (jllgner fc Nelson 
20061 ) and our model could - possibly with an ad-hoc re- 
duction factor for the drift speed - provide new insight in 
this kind of modeling. 

Our disk model involves a constant threshold fragmen- 
tation velocity which can be put into question. Laboratory 
experiments show that this threshold velo c ity is dependent 
on p article size ( Blum fc MuenchI Il993t IParaskov et al.l 
20071 ) and we indeed work on that topic to investigate 
if more realistic aggregate collision models predict a dif- 
ferent disk evolution (Brauer et al. in prep.). 

Finally, current state-of-the-art models of disk struc- 
ture and their spectra and images relies on ad-hoc pre- 
scriptions of the dust spatial and size distribution. Models 
of the kind described here will be linked to radiative trans- 
fer calculations to investigate the effect of grain evolution 
on disk structure in the near future. 
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Appendix A: Coagulation algorithms 

A.l. Podolak algorithm 

In this section we present an algorithm that was first used be lKovetz fc Olundl ( 1969 ) in meteorological science. This 



algorithm conserves mass and particle number density. It captivates by its simplicity and it is comparatively easy to 
implement into a computer code. 

Let us assume a mass grid rrii and a given spatial number density A'i of the particles with mass mi. Two particles 
of mass TTii and ttij coagulate with a coagulation rate Qij, i.e. the number of coagulation events per time, which is 
given by 

Qij = iViTVji^ij . (A.l) 

The quantity Kij denotes the coagulation kernel of the particles in the mass bins i and j. It is given by the product of 
the collisional cross section Uij between the particles and their relative velocity wy . High values of this kernel correspond 
to high collision rates which means that the particles are subject to a fast growth process. Low values of this quantity 
imply slow growth rates. 

Now, an important issue appears if non-linear mass grids are considered. The resulting mass of the coagulation, 
i.e. the mass m ^ mi + rrij, does not neccessarily match with any of the mass grid points. This means that in general 
no mass grid point rus can be found which satisfies rus = mi + rrij. Therefore, we have to divide the coagulating mass 
between the nearest mass grid points in some sensible way. 

We assume that the nearest neighbours are given by Trim < m < nin- With a linear ansatz we split the coagulation 
rate Qij into a coagulation rate for the mass m^ and a coagulation rate for the mass rrin, 

Qm = e Qij and Q„ = (1 - e) Qij . (A.2) 

The number density is a conserved quantity in this algorithm since Q^ + Qn ~ Qij . Much more important than the 
conservation of number density is the conservation of mass. We can enforce this fundamental conservation principle 
by setting 

Qm"lm + Qn'Tln = Qij (mj + TOj) . (A. 3) 

The last expression defines the value of e which was a free parameter till now. Inserting Eqs. (|A.2p into Eq. (jA.Sp we 

find that e can be written as 

_ mn — (m\ + mj ) 



(A.4) 



''n ^'^m 



One may translate these e parameters for every coagulation process between particles of species i and j into certain 
coefficients Cijk so that the coagulation equation can be expressed in the form 



^k = ^E^y^yk-J^Qik. (A.5) 

The quantity C is then given by 



2 



e if my is the largest mass grid point < vn,\ + toj 

Cijk — ^ 1 ^ ^ if '^k is the smallest mass grid point > m\ + m-^ 

otherwise 

In general, more than 90% of the elements of the matrix C are zero. Therefore, a lot of computer calculation time can 
be saved if only the non-zero elements in the last expression are summed up. 

A.l. Modified Podolal< algorittim 

In simulations for protoplanetary disks, we are interested in particle growth from sub-micrometer in size, i.e. particles 
that are part of the interstellar medium, up to several hundred meters. This means that particle masses ranging from 
10~i2 g to 10^^ g are considered, which corresponds to 24 orders of magnitude in mass. 

Now, a major problem appears. Intrinsic computer variables, so-called double precision variables, have an accuracy 
only up to 14 digits. Since we are interested in particle growth by more than 20 orders of magnitude, the accuracy needed 
for coagulation simulations exceeds the accuracy provided by the computer. One solution could be the introduction of 
quadrupole precision computer variables but this would slow down the simulation speed significantly. This accuracy 
issue leads to certain problems for example the violation of mass conservation or simulation crashes. In order to perform 
the coagulation simulations anyway we have to analytically reformulate the Podolak algorithm at certain points of the 
numerical scheme. 
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For this purpose we first have to introduce a number Cg in the following way. We consider the neighboring mass 
grid points mk-i and rrik- The number Ce is then defined in a way that the inequality 

mk-i + rrii < my^ (A. 6) 

holds for any i which satisfies the condition i < k — Ce- In general, the value of Ce is dependent on the index k. 

Now, we reformulate the algorithm. We formally separate the diagonal elements of Eq. (jA.SP from the non-diagonal 
elements. 

k k i-1 N 

^k = 2 S ^^'^iki^ii + E E ^i^J^iJk^iJ - E ^J^ki^Jk • (A.7) 

i=l i=l j=l j=l 

In the next step, we consider the second and third summand of the last expression. In the second term we separate 
the case i = k which leads to 

k-l i-1 k+l-Ce N 

EE^i^jC'ijki^iJ+ E ^kA^jCkjk/^kj-E^J^ki^Jk- (A-8) 

i=l j = l j = l j 

This can be rewritten as 

k-l i-1 k+l-Ce N 

EE^i^J^ykifu+ E (^kA^jCkjkA\j-7Vj7Vki^jk)- E ^J^kifjk 

i=l j = l j = l j=k+2-c, 

k-l i-1 k+l-Ce N 

= EE^i^jC'ijki^y- E ^kA^j-^kj ""i - E ^j^ki^ik 

■ 1-1 -1 '"-k+l '"-k ■ 1 I o 

1=1 j = l j = l j=k+2-Ce 

k-l i-1 N 

= E E ^i^jC-ijki^ij + E ^kiVji^kji^jk , (A.9) 

i=l j=l j=l 

where the matrix D is given by 

if J < k + 1 — Ce and 



' I -1 ifj>k+l-Ce. 

The new coagulation equation now reads 

k k-l i-1 AT 

^k - 2 S ^i" ^"l^^" + E E ^'^j ^ijki^ij + E ^kA^j i^kj i^jk . (A. 10) 

i=l i=l j=l j=l 

This was one part of rewriting the algorithm. For the other part we regard the second term of the last expression, 
especially the term i = k — 1. We can rewrite this term as follows, 

k-2 

E^k-iAjCk-ij,kifk-ij 

k-Cc fc-2 

= E ^k-iAjCk-i,j,ki^k-i.j + E ^k-iAfjCk-ij^kATk-ij 

j = l k-Ce + l 

k-Ce k-2 



E^k-iA'j ^ ^l j^k-i.j+ E A^k-iA'jCk-i,j.kifk-ij 

J = l k-Ce + 1 

k-2 

E^k-iA^ji^k-iji^jk 

j=l 

N N / 3\ 

E E NiN^K.^E^s+iO k - j - 2 J -^i.k-i . (A.ll) 

i=i i=i ^ ^ 



In this Equation the matrix E is given by 

if j < k — Ce and 
' (?7ik+i -m-i - mk-i) if j > k - Ce . 



J mk-nik-l 
-C-ik — 



^jk 



, __ mj+nik-l- "Ik 
mk+i-nik 
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With these two reformulations the coagulation equation can be written in the form 



(A.12) 



where the final coagulation matrix M is given by 

Mj = \hCiik + Cijke (k - i - ^ j (i - J - + '^iki^ji + '5i,k-i£;j,i+ie \^-i-\ 

In this expression Q(x) denotes the Heaviside distribution, which is zero for a; < and unity for x > 1. 



Appendix B: Vertical integration 

Coagulation and fragmentation are local processes. This means that the equations described in the last section have to 
be solved at every point in space. The more space grid points are considered the more time-consuming the computer 
simulations become. However, under certain conditions the situation simplifies. In the following we will describe a 
scheme that can save a remarkable amount of computational time. 

We consider the coagulation (fragmentation) equation at a certain space point z-p 



7Vk(^p) = ^Gijk(zp)M(zp)iVj(zp) 



(B.l) 



Since we are interested in particle growth in protostellar disks we can adapt the number densities to this special 
problem. We assume that at any given time the vertical particle distribution of any given particle size is given by a 
settling-mixing equilibrium distribution (cf. Eg. llSp . This leads to a density N^ of a particle of size Oi which depends 
on the height above the midplane z as 



N-Az) 



2TThi 



exp 



1 ± 

2 V/ii 



(B.2) 



In this expression the variable hi denotes the dust scale height of the particles with mass nii. The quantity LUi is the 
surface number density of the particles with that certain mass. Inserting Eq. (jB.2p into Eq. (jB.ip and integrating over 
height above the midplane z yields 






hi 



exp 



A^ 



If we define 



Gijk = Yl 



Gij]i(Zp) 
2T:hihj 



exp 



exp 



^ 



Az 



the integrated coagulation equation can be written as 



Wk 



ij 



ujiUjjGijk 



(B.3) 



(B.4) 



(B.5) 



In this way we have integrated the z-dimension out without a single approximation only with the assumption that the 
vertical redistribution goes faster than the coagulation/fragmentation. This reformulation of the coagulation equation 
has an obvious advantage, namely instead of solving the coagulation equation at every point in z the last expression 
enables us to solve the equation for every height above the midplane at the same time. If we assume a vertical grid 
with 60 grid points the vertical integration speeds up the computer simulation routine by a factor of 60. 



Appendix C: Implicit differencing 

If fragmentation is included in the simulations then the limiting time step for the coagulation/fragmentation process 
tends to be small. Fragmentation leads to a permanent amount of small particles. Small particles, however, are 
associated with short time scales. Taking these short time scales into account, the time step of the numerical simulation 
can not be chosen to be very large. This argumentation only holds for explicit numerical solvers. For this reason we have 
implemented an implicit solver for the coagulation/fragmentation equation which we will describe in the following. 
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The coagulation/fragmentation equation can be written in the form 
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(C.l) 



where / denotes the particle distribution vector on the mass grid and the function F describes the time evohition. In 
one time step At at a certain time t, we now want to calculate the new particle distribution fn — f{t + At) from the 
old distribution fo = J{t). Therefore we rewrite Eq. (|C.1|1 as 



e^AtFiU), 



(C.2) 



where t = fn — fo and fi = ^/o + (1 — C)fn- The time evolution of the function / with ^ = 1 is called explicit, while the 
time evolution with ^ = is usually called implicit. Choosing ^ = in our case, we can perform a Taylor expansion of 
the right-hand side of Eq. (jC.2|l which leads to 



e = AtF{fo) + AtJl . 
The Matrix J denotes the Jacobi matrix which is definded as Jy = dFi/dfj. Solving Eq. (jC.Sp for e leads to 



l-AtJ 



AtFifo) 



(C.3) 



(C.4) 



Hence, the evolution of the implicit time step reduces to a matrix inversion which can be done easily. 
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Fig. D.l. Test for the radial drift routine. Coagulation is 
neglected in this simulation. The lines denote the parti- 
cle position in the disk after 10^ yrs of radial drift from 
the initial radial position ro = 5.53 AU for different dust 
particle sizes. 



Appendix D: Tests 

This part of the appendix considers a comparison between 
the model described in Section [5] and a model which in- 
volves the following approximation. The particle distribu- 
tions which can be seen for example in Fig. [H] are fairly 
narrow peaks. Hence, it is suggestive to approximate the 
particle distribution by a monodisperse distribution, i.e. a 
distribution with only one single particle size. In this case 
the coagulation equation simplifies enormously. We will 
compare these two different models by considering coagu- 
lation due to Brownian motion and turbulent coagulation. 
Also the routine for radial motion will be checked against 
this simplified model. 

D.l. Simplified model and numerical setup 

We assume that a certain amount of equally-sized dust 
particles with radius a is located at a single radius r in 
the disk. The surface density of these particle is given by 
Ed • The dust particles are vertically distributed according 
to Eq. (fTSj) . For the test cases we ignore the radial motion 
of the gas and the radial turbulent diffusion of the dust 
so that the radial motion of the dust is given by Eq. PT|) . 
The assumption of a monodisperse distrib ution leads to 
a coa gulation equation given by Eq. ([34)) (JKornet et al.l 
l200lh . In the test cases we will use a surface density of 
the dust Ed = 5.26 x 10~^ g/cm^ and an initial location 
in the disk given by ro = 5.53 AU. Every other parameter 
was already mentioned in Section [2 



D.2. Radial drift 

First, we check the radial drift routine. Any coagulation 
is neglected. With the values mentioned above we let the 



Fig. D.2. Test for the coagulation routine. This plot 
shows the particle distribution in the full model and the 
monodisperse model for coagulation due to Brownian mo- 
tion at 4 different times. The stars '•' denote the particle 
size in the monodisperse model. The largest discrepancy 
in particle size between the two models is a factor of 1.6 
in particle size a. 



dust particles drift for 10^ yrs and plot their radial po- 
sition in the disk for different particle sizes. The results 
of this simulation are shown in Fig. ID. II In this figure 
the solid line correponds to the dominant particle size of 
the dust distribution in the full model. The dashed line 
denotes the result of the monodisperse model. Any dis- 
continuous effects are due to the radial grid. 

D.3. Coagulation 

Now, we consider the coagulation of the dust while the 
radial motion of the dust is neglected. We investigate the 
dust particle coagulation at rg = 5.53 AU in the disk and 
we first focus on coagulation due to Brownian motion. The 
results of both models are shown in Fig. ID. 21 The stars 
't^t' denote the particle size in the monodisperse model. 
The largest discrepancy in particle size between the two 
models is a factor of ^ 1.6 after 10^ yrs. 

The results of the same simulation but now with par- 
ticle growth due to turbulent coagulation included are 
shown in Fig. ID. 31 The largest discrepancy in this case in 
particle size between the two models is a factor of ~ 2.7 
in radius for t = W^ yrs. 

However, if the Stokes number of the dust particles is 
smaller than unity, then the particle size of both mod- 
els after a certain time can differ in more than one order 
of magnitude. This is due to the following reason. If the 
Stokes number is smaller than unity then the relative tur- 
bulent velocity of the dust follows Av ex \/St- With the 
coagulation equation ([34]) this leads to a oc a and its solu- 
tion a ex exp(t). In the particular case of turbulent coag- 
ulation and St < 1, the particle size as a function of time 
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Fig. D.3. Test for the coagulation routine. This plot 
shows the particle distribution in the full model and the 
monodisperse model for coagulation due to Brownian mo- 
tion and turbulence in the disk at 3 different times. The 
stars '•' denote the particle size in the monodisperse 
model. The largest discrepancy in particle size between 
the results that are shown in this figure is a factor of ^ 2.7. 



is not given by a polynomial expression, i.e. an expression 
of the form a (xV , but an exponential function. If a time 
evolution follows a polynomial law then two different but 
similar initial conditions will lead to different but simi- 
lar particle sizes after any time. Let us consider the time 
evolution 



a{t) = 



ct 



1/7 



(D.l) 



First, the initial particle size oo gets unimportant for 
large t. Second, two different c-parameters, i.e. ci and C2, 
will lead to particle sizes which will differ by a factor of 
(ci/c2)'*^. This factor is constant and does not increase. 
On the other hand, if the particles grow exponentially, 



a{t) = oo exp(ci), 



(D.2) 



then the initial particle size ao will always play a role. 
Moreover, two different c- values, i.e. Ci and C2, will lead to 
particle sizes which will differ by a factor of exp(cit — C2t) . 
This factor increases in time. For example, in the case of 
turbulent coagulation (cf. eq. [32]) the parameter c is given 
by eo^^k- An initial dust-to-gas ratio of 1% leads to a = 
ape « 2.72ao after 100 orbits. An initial dust-to-gas ratio 
of 2% after 200 orbits already imphes a = aoe* « 54.6ao 
which is a factor of 20 larger. In the case of turbulent 
coagulation and St smaller than unity, small changes in 
the intial conditions lead to large differences in the growth 
behaviour. Hence, a change from a monodisperse model to 
a model with a whole particle dispersion presumably leads 
to similar effects. 



Appendix E: Collision rates 

In this part of the Appendix, we investigate the equilib- 
rium particle distribution between coagulation and frag- 
mentation after 10'^ yrs of disk evolution at 1 AU in 
the disk. In this calculation the fragmentation velocity is 
Vf — l{fi cm/s and the fragmentation parameter ^ is 1.83. 
We adopt a disk mass of 10^^ A/^, a turbulent a- value of 
10"'^ and an initial dust-to-gas ratio of 10^^. The crater- 
ing parameter x = 0.5 and we adopt ■0 — 2. We focus on 
the collision rates between particles of different sizes and 
how important these collisions are for particle growth. 

We define the collision rates -R(ri, r2) as the number of 
collisions per second between particles of radius ri and r2 
in a vertical column with a cross section of 1 cm^. These 
collision rates are shown in Fig. IE. II The collision rates for 
particles of equal size decrease dramatically with increas- 
ing particle radius. The collision rates of micrometer-sized 
particles and the collision rates of mm-sized particles differ 
by more than 10 orders of magnitude. While the collision 
rate for 1 /im particles is ^ lO''"^ cm~^yrs~^, the collision 
rate for 10 /xm particles is already more than two orders 
of magnitude lower. The smallest particles have the high- 
est collision rates. Since the coagulation probability for 
fjLTn particles is rather high, this leads to relatively short 
coagulation time scales for small particles. The collision 
rate is 1 cm^^yrs^^ for 10^ /im particles. For particles, 
which are one order of magnitude in radius larger, this rate 
has already dropped to 10"^ cm^^yrs^^. However, the be- 
haviour of the collision rates between non-equal sizes par- 
ticles is different. For example, the collision rate between 
'^ /^m-sized dust particles and larger particles does hardly 
change over a wide range. 

The collision rates, which are shown in Fig. IE. II do 
not provide information about the importance of collisions 
for particle growth. Not only the number of collisions per 
time is important for particle coagulation, but also the 
mass of the particles itself. For this reason, we calculate 
the relative mass gain of a single particle with radius ri . If 
the number of all particles of radius ri in a vertical column 
is given by A^i, then in average every single particle of 
size ri collides with R{ri,r2)/Ni particles of size r2 per 
second. We assume that the sticking probability is unity. 
This means that the mass, which the particle of radius ri 
sweeps up in time, is given by R{ri, r-2)m2lN\. Its relative 
mass gain rate is then given by 



2ri(?-2) 



R{ri,r2) 7712 



A^i 



77^1 



(E.l) 



This quantity is shown in Fig. IE. 21 Since larger particles 
sweep up smaller particles the mass gain rates S are not 
shown for particle sizes r2 > ri. 

We find that the collisions which are most impor- 
tant for particle growth are collisions between equal-sized 
particles. This statement holds for particles smaller than 
~ 0.5 mm in size. If the particle with radius ri is larger 
than this value then collisions with r2 ~ 0.5 mm sized 
particles are most important for the growth of the dust. 
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Fig. E.l. This plot shows the verticaUy integrated num- 
ber of coUisions per time at 1 AU in the disk as discussed 
in Sec. [El This calculation is based on the particle distri- 
bution after 10'^ yrs of disk evolution shown in Fig. [TS] In 
this simulation we adopted a disk mass of 10~^ A/^ and a 
turbulent a parameter of 10~^. 
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Fig. E.2. This plot shows the relative mass gain per year 
at 1 AU in the disk. This calculation is based on the par- 
ticle distribution after 10^ yrs of disk evolution. The disk 
mass is 10~^ Mi, and the turbulent a parameter is 10"'^. 
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Variable Explanation 



a 


radius of the particle 


m 


mass of the particle 


ps 


solid material density of the particle 



distance to the central star from a point in the midplane 

height above the midplane 

temperature 

isothermal soundspeed 

Kepler frequency, Kepler velocity 

gas scale height 

dust scale height 

surface density of the gas and the dust 

gas and dust density 

initial and current dust-to-gas ratio 

inner radius of the disk 

outer radius of the disk 

mass of the central star 

mass of the disk 

tin maximum radial drift velocity 

Q, q turbulence parameters 

St Stokes number of the particle 

Dg , _Dd diffusion coefficients of gas and dust 

iidust radial dust drift velocity due to gas drag 

iigas radial gas accretion velocity 

iiju^st total radial dust velocity 

^ slope of the particle distribution after fragmentation 

X relative amout of mass removed from the target particle by cratering 

Vf threshold fragmentation velocity 

Table E.l. Important variables used in the course of this paper. 
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